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Abstract 


This paper introduces time-continuous numerical schemes to simulate sto¬ 
chastic differential equations (SDEs) arising in mathematical finance, population 
dynamics, chemical kinetics, epidemiology, biophysics, and polymeric fluids. These 
schemes are obtained by spatially discretizing the Kolmogorov equation associated 
with the SDE in such a way that the resulting semi-discrete equation generates a 
Markov jump process that can be realized exactly using a Monte Carlo method. 
In this construction the spatial increment of the approximation can be bounded 
uniformly in space, which guarantees that the schemes are numerically stable for 
both finite and long time simulation of SDEs. By directly analyzing the generator 
of the approximation, we prove that the approximation has a sharp stochastic Lya¬ 
punov function when applied to an SDE with a drift field that is locally Lipschitz 
continuous and weakly dissipative. We use this stochastic Lyapunov function to 
extend a local semimartingale representation of the approximation. This extension 
permits to analyze the complexity of the approximation. Using the theory of semi¬ 
groups of linear operators on Banach spaces, we show that the approximation is 
(weakly) accurate in representing finite and infinite-time statistics, with an order 
of accuracy identical to that of its generator. The proofs are carried out in the con¬ 
text of both fixed and variable spatial step sizes. Theoretical and numerical studies 
confirm these statements, and provide evidence that these schemes have several ad¬ 
vantages over standard methods based on time-discretization. In particular, they 
are accurate, eliminate nonphysical moves in simulating SDEs with boundaries (or 
confined domains), prevent exploding trajectories from occurring when simulating 
stiff SDEs, and solve first exit problems without time-interpolation errors. 
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CHAPTER 1 


Introduction 

1.1. Motivation 

Stochastic differential equations (SDEs) are commonly used to model the effects 
of random fluctuations in applications such as molecular dynamics HM im n 
134] , molecular motors [6], microfluidics [8311129) . atmosphere/ocean sciences [95j . 
epidemiology iMi mu, population dynamics [3 m], and mathematical finance 
[Mlilll]. In most of these applications, the SDEs cannot be solved analytically 
and numerical methods are used to approximate their solutions. By and large the 
development of such methods has paralleled what has been done in the context of 
ordinary differential equations (ODEs) and this work has led to popular schemes 
such as Euler-Maruyama, Milstein’s, etc |72l \Tm [TOT] . While the design and 
analysis of these schemes differ for ODEs and SDEs due to peculiarities of stochastic 
calculus, the strategy is the same in both cases and relies on discretizing in time 
the solution of the ODE or the SDE. For ODEs, this is natural: their solution is 
unique and smooth, and it makes sense to interpolate their path at a given order of 
accuracy between snapshots taken at discrete times. This objective is what ODE 
integrators are designed for. The situation is somewhat different for SDEs, however. 
While their solutions (in a pathwise sense) may seem like an ODE solution when 
the Wiener process driving them is given in advance, these solutions are continuous 
but not differentiable in general, and an ensemble of trajectories originates from 
any initial point once one considers different realizations of the Wiener process. 

The simulation of SDEs rather than ODEs presents additional challenges for 
standard integration schemes based on time-discretization. One is related to the 
long time stability of schemes, and their capability to capture the invariant distri¬ 
bution of the SDE when it exists. Such questions are usually not addressed in the 
context of ODEs: indeed, except in the simplest cases, it is extremely difficult to 
make precise statements about the long time behavior of their solutions since this 
typically involves addressing very hard questions from dynamical systems theory. 
The situation is different with SDEs: it is typically simpler to prove their ergodicity 
with respect to some invariant distribution, and one of the goals of the integrator is 
often to sample this distribution when it is not available in closed analytical form. 
This aim, however, is one that is difficult to achieve with standard integrators, be¬ 
cause they typically fail to be ergodic even if the underlying SDE is. Indeed the 
numerical drift in standard schemes can become destabilizing if the drift field is 
only locally Lipschitz continuous, which is often the case in applications. These 
destabilizing effects typically also affect the accuracy of the scheme on finite time 
intervals, making the problem even more severe. Another important difference be¬ 
tween SDEs and ODEs is that the solutions of the former can be confined to a 
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certain region of their state space, and it is sometimes necessary to impose bound¬ 
ary conditions at the edge of this domain. Imposing these boundary conditions, or 
even guaranteeing that the numerical solution of the scheme remains in the correct 
domain, is typically difficult with standard integration schemes. 

In view of these difficulties, it is natural to ask whether one can take a differ¬ 
ent viewpoint to simulating SDEs, and adopt an approach that is more tailored to 
the probabilistic structure of their solution to design stable and (weakly) accurate 
schemes for their simulation that are both provably ergodic when the underlying 
SDE is and faithful to the geometry of their domain of integration. The aim of 
this paper is to propose one such novel strategy. The basic idea is to discretize 
their solution in space rather than in time via discretization of their infinitesi¬ 
mal generator. Numerous approximations of this second-order partial differential 
operator are permissible including finite difference or finite volume discretizations 
|3T1 fM fT03l fl04l [85]. As long as this discretization satisfies a realizability 
condition, namely that the discretized operator be the generator of a Markov jump 
process on a discrete state space, it defines a process which can be exactly simu¬ 
lated using the Stochastic Simulation Algorithm (SSA) [371131L I142L 185] . Note 
that this construction alleviates the curse of dimensionality that limits numerical 
PDE approaches to low dimensional systems while at the same time permitting to 
borrow design and analysis tools that have been developed in the numerical PDE 
context. The new method gets around the issues of standard integrators by permit¬ 
ting the displacement of the approximation to be bounded uniformly and adaptively 
in space. This feature leads to simple schemes that can provably sample the sta¬ 
tionary distribution of general SDEs, and that remain in the domain of definition 
of the SDE, by construction. In the next section, we describe the main results of 
the paper and discuss in detail the applications where these SDE problems arise. 


1.2. Main Results 

If the state space of the approximating Markov jump process, or approxima¬ 
tion for short, is finite-dimensional, then a large collection of analysis tools can be 
transported from the theory of numerical PDE methods to assess the stability and 
accuracy of the approximation. However, if the state space is not finite-dimensional, 
then the theoretical framework presented in this paper and described below can be 
used instead. It is important to underscore that the Monte Carlo method we use to 
realize this approximation is local in character, and hence, our approach applies to 
SDE problems which are unreachable by standard methods for numerically solving 
PDEs. In fact, the approximation does not have to be tailored to a grid. Leveraging 
this flexibility, in this paper we derive realizable discretizations for SDE problems 
with general domains using simple finite difference methods such as upwinded and 
central difference schemes [mllMl l42]. though we stress that other types of dis¬ 
cretizations that satisfy the realizablility condition also fit our framework. 

Since an outcome of this construction is the generator of the approximation, 
the proposed methods have a transparent probabilistic structure that make them 
straightforward to analyze using probabilistic tools. We carry out this analysis in 
the context of an SDE with unbounded coefficients and an unbounded domain. 
Our structural assumptions permit the drift field of the SDE to be locally Lipschitz 
continuous, but require the drift field to satisfy weak dissipativity and polynomial 
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growth conditions. Our definition of stability is that the approximation has a sto¬ 
chastic Lyapunov function whenever the underlying SDE has one. This definition is 
natural when one deals with SDEs with locally Lipschitz drift fields. Indeed, global 
existence and uniqueness theorems for such SDEs typically rely on the existence of 
a stochastic Lyapunov function [67ll97|. 

In this context, here are the main results of the paper regarding realizable 
discretizations with gridded state spaces. We emphasize that we do not assume 
that the set of all grid points is finite. 

(1) We use Harris Theorem to prove that the approximation is geometrically 
ergodic with respect to an invariant probability measure [1061 m 116]. 
To invoke Harris theorem, we prove that the approximation preserves a 
sharp stochastic Lyapunov function from the true dynamics. Sharp, here, 
means that this stochastic Lyapunov function dominates the exponential 
of the magnitude of the drift field, and can therefore be used to prove that 
a large class of observables (including all moments of the approximation) 
are integrable. 

(2) In addition to playing an important role in proving existence of an in¬ 
variant probability measure and geometric ergodicity, this stochastic Lya¬ 
punov function also helps solve a martingale problem associated to the 
approximation [Il[TT6l[32]. This solution justifies a global semimartin¬ 
gale representation of the approximation, and implies that Dynkins for¬ 
mula holds for the conditional expectation of observables that satisfy a 
mild growth condition. We apply this semimartingale representation to 
analyze the complexity of the approximation. Specifically, we obtain an 
estimate for the average number of computational steps on a finite-time 
interval in terms of the spatial step size parameter. 

(3) We use a variation of constants formula to quantify the global error of 
the approximation ffni This formula is a continuous-time analog of the 
Talay-Tubaro expansion of the global error of a discrete-time SDE method 
[HI]. An analysis of this formula reveals that the order of accuracy of 
the approximation in representing finite-time conditional expectations and 
equilibrium expectations is given by the order of accuracy of its generator. 
Geometric ergodicity and finite-time accuracy imply that the approxima¬ 
tion can accurately represent long-time dynamics, such as equilibrium 
correlation functions [14j . 

This paper also expands on these results by considering space-discrete generators 
that use variable spatial step sizes, and whose state space may not be confined to a 
grid. In these cases, we embed the ‘gridless’ state space of the approximation into 
K", and the main issues are: 

• the semigroup of the approximation may not be irreducible with respect 
to the standard topology on K"; and 

• the semigroup of the approximation may not even satisfy a Teller property 
because its associated generator is an unbounded operator and the sample 
paths of its associated process are discontinuous in space ungdin]. 

To deal with this lack of irreducibility and regularity, we use a weight function to 
mollify the spatial discretization so that the resulting generator is a bounded linear 
operator. This mollification needs to be done carefully in order for the approxima¬ 
tion to preserve a sharp stochastic Lyapunov function from the true dynamics. In 
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short, if the effect of the mollification is too strong, then the approximation does 
not have a sharp stochastic Lyapunov function; however, if the mollification is too 
weak, then the operator associated to the approximation is unbounded. 

The semigroup associated to this mollified generator is Feller, but not strongly 
Feller, because the driving noise is still discrete and does not have a regularizing 
effect. However, since the process is Feller and has a stochastic Lyapunov function, 
we are able to invoke the Krylov-Bogolyubov Theorem to obtain existence of an 
invariant probability measure for the approximation ffnj. This existence of an 
invariant probability measure is sufficient to prove accuracy of the approximation 
with respect to equilibrium expectations. To summarize, here are the main results 
of the paper regarding realizable discretizations with gridless state spaces. 

(1) The approximation is stable in the sense that it has a sharp stochastic 
Lyapunov function, and accurate in the sense that the global error of the 
approximation is identical to the accuracy of its generator. The proof of 
both of these properties entails analyzing the action of the generator on 
certain test functions, which is straightforward to do. 

(2) A semimartingale representation for suitable observables of the approx¬ 
imation is proven to hold globally. This representation implies Dynkins 
formula holds for this class of observables. 

(3) The approximation has an invariant probability measure, which may not 
be unique. This existence is sufficient to prove that equilibrium expecta¬ 
tions of the SDE solution are accurately represented by the approximation. 

To be clear, the paper does not prove or disprove uniqueness of an invariant prob¬ 
ability measure of the approximation on a gridless state space. One way to resolve 
this open question is to show that the semigroup associated to the approximation 
is asymptotically strong Feller |44) . 

Besides these theoretical results, below we also apply the new approach to a 
variety of examples to show that the new approach is not only practical, but also 
permits to alleviate several issues that commonly arise in applications and impede 
existing schemes. 

(1) Stiffness problems. Stiffness in SDEs can cause the numerical trajectory to 
artificially “explode.” For a precise statement and proof of this divergence 
in the context of forward Euler applied to a general SDE see |59j . This 
numerical instability is a well-known issue with explicit discretizations 
of nonlinear SDEs [TM ism no8i na [60], and it is especially acute if 
the SDE is stiff e.g. due to coefficients with limited regularity, which is 
typically the case in applications. Intuitively speaking, this instability is 
due to explicit integrators being only conditionally stable, i.e. for any fixed 
time-step they are numerically stable if and only if the numerical solution 
remains within a sufficiently large compact set. Since noise can drive 
the numerical solution outside of this compact set, exploding trajectories 
occur almost surely. The proposed approximation can be constructed 
to be unconditionally stable, since its spatial step size can be bounded 
uniformly in space. As a result, the mean time lag between successive 
jumps adapts to the stiffness of the drift coefficients. 

(2) Long-time issues. SDEs are often used to sample their stationary dis¬ 
tribution, which is generically not known, meaning there is no explicit 
formula for its associated density even up to a normalization constant. 
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Examples include polymer dynamics in extensional, shear, and mixed 
planar flows m HH HH Ull |126]; polymers in confined flows 

[I9l [631 [641 [651 fl43l 11441 [501 [491 [481 fSTlIin 1^1146] : polymer confor¬ 
mational transitions in electrophoretic flows IM EH IM SMI HI m] ; 
polymer translocation through a nanopore |91L I88L IIOIL I56L |4]; and 
molecular systems with multiple heat baths at different temperatures 
[301 \29[ [IT71 [991 [451 [871 [T] as in temperature accelerated molecular 
dynamics |99L I139L [1]. Unknown stationary distributions also arise in 
a variety of diffusions with nongradient drift fields like those used in at¬ 
mosphere/ocean dynamics or epidemiology. The stationary distribution 
in these situations is implicitly defined as the steady-state solution to a 
PDE, which cannot be solved in practice. Long-time simulations are there¬ 
fore needed to sample from the stationary distribution of these SDEs to 
compute quantities such as the equilibrium concentration of a species in 
population dynamics or, the polymer contribution to the stress tensor 
or equilibrium radius of gyration in Brownian dynamics. This however 
requires to use numerical schemes that are ergodic when the underlying 
SDEs they approximate are, which is by no means guaranteed. We show 
that the proposed approximation is long-time stable as a consequence of 
preserving a stochastic Lyapunov function from the true dynamics. In 
fact, this stochastic Lyapunov function can be used to prove that the ap¬ 
proximation on a gridded state space is geometrically ergodic with respect 
to a unique invariant probability measure nearby the stationary distribu¬ 
tion of the SDE, as we do in this paper. 

(3) Influence of boundaries. SDEs with boundaries arise in chemical kinet¬ 

ics, population dynamics, epidemiology, and mathematical finance, where 
the solution (e.g., population densities or interest rates) must always be 
non-negative [#l HH [1221 [221 [2l [TSl [gm [5]. Boundaries also arise 
in Brownian dynamics simulations. For example, in bead-spring chain 
models upper bounds are imposed on bond lengths in order to model the 
experimentally observed extension of semi-flexible polymers In this 

context standard integrators may produce updated states that are not 
within the domain of the SDE problem, which in this paper we refer to 
as nonphysical moves. The proposed approximation can eliminate such 
nonphysical moves by allowing the spatial step size to be set such that no 
update lies outside the domain of definition of the SDE. 

(4) First Exit Problems. It is known that standard integrators do not work 
well in simulating first exits from a given region because, simply put, they 
fail to account for the probability that the system exits the region in be¬ 
tween each discrete time step [3811391140) . Building this capability into 
time integrators is nontrivial, e.g., one technique requires locally solving 
a boundary value problem per step, which is prohibitive to do in many 
dimensions [96j . By using a grid whose boundary conforms to the bound¬ 
ary of the first exit problem, the proposed approximation does not make 
such time interpolation errors. 

In addition, the examples illustrate other interesting properties of realizable dis¬ 
cretizations such as their accuracy with respect to the spectrum of the infinitesimal 
generator of non-symmetric diffusions. 
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1.3. Relation to Other Works 

When the diffusion matrix is diagonally dominant, H. Kushner developed real¬ 
izable finite difference approximations to SDEs and used probabilistic methods to 
prove that these finite difference schemes are convergent f76l[80ll8nPf71Pf8lP79l . 
When the SDE is self-adjoint and the noise is additive, realizable discretizations 
combined with SSA have also been developed, see [ 3 n[T 4 ^[T 03 llT 04 l[ 85 ] . and to 
be specific, this combined approach seems to go back to at least m Among these 
papers, the most general and geometrically flexible realizable discretizations seem 
to be the finite volume methods developed in |85j . This paper generalizes these 
realizable discretizations to diffusion processes with multiplicative noise. 

Our approach is related to multiscale methods for approximating SDEs that 
leverage the spectrum of suitably defined stochastic matrices mg [Till US usi 
[TM [IMI [TM [IM m- However, to be clear, this paper is primarily about 
approximating the SDE solution, not coarse-graining the dynamics of the SDE. 

Our approach also seems to be foreshadowed by Metropolis integrators for 
SDEs, in the sense that the proportion of time a Metropolis integrator spends in a 
given region of state space is set by the equilibrium distribution of the SDE, not the 
time step size mni mi [n mini aa m ■ The idea behind Metropolis integrators 
is to combine an explicit time integrator for the SDE with a Monte Carlo method 
so that the resulting scheme is ergodic with respect to the stationary distribution 
of the SDE and finite-time accurate |13j . However, a drawback of these methods is 
that they are typically limited to diffusions that are self-adjoint m, whereas the 
methods presented in this paper do not require this structure. Moreover, Metropolis 
integrators are typically ergodic with respect to the stationary distribution of the 
SDE, but they may not be geometrically ergodic when the SDE is [TlOllTllfTT]. 
Intuitively speaking, the reason Metropolis integrators lack geometric ergodicity 
is that they use as proposal moves conditionally stable integrators for the SDE. 
Hence, outside the region of stability of these integrators the Metropolis rejection 
rate tends to one. In these unstable regions. Metropolis integrators may get stuck, 
and as a consequence, typically do not have a stochastic Lyapunov function. In 
contrast, the proposed approximations do have a stochastic Lyapunov function for 
a large class of diffusion processes. 

Our approach is also related to, but different from, variable time step integra¬ 
tors for SDEs [36111321182| . These integrators are not all adapted to the filtration 
given by the noise driving the approximation. Let us comment on the adapted 
(or non-anticipating) technique proposed in [82] , which is based on a simple and 
elegant design principle. The technique has the nice feature that it does not require 
iterated stochastic integrals. The method adjusts its time step size at each step 
according to an estimate of the ‘deterministic accuracy’ of the integrator, that is, 
how well the scheme captures the dynamics in the zero noise limit. The time step 
size is reduced until this estimate is below a certain tolerance. The aim of this 
method is to achieve stability by being deterministically accurate. In fact, as a 
consequence of this accuracy, the approximation ‘inherits’ a stochastic Lyapunov 
function from the SDE and is provably geometrically ergodic. This paper [82j also 
addresses what happens if the noise is large compared to the drift. The difference 
between these variable time step integrators and the proposed approximation is 
that the latter sets its mean time step according to the Kolmogorov equation and 
achieves stability by bounding its spatial step size. 
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1.4. Organization of Paper 

The remainder of this paper is structured as follows. 


Chapter details our strategy to spatially discretize an SDE problem. 

defines what it means for a spatial discretization to be realizable, 
describes the differences and pros/cons between a gridded and gridless 


S2.1 


? 2.2 


state space. 


^2.3 derives realizable discretizations in one dimension using an upwinded 
finite difference, finite volume, and a probabilistic discretization. 

§2.4| derives realizable discretizations in two dimension using an upwinded 
finite difference method. 


^2.5 presents realizable discretizations in n-dimension, which are based on 


central and upwinded finite differences on a gridless state space. 

§2.6| discusses and analyzes the scaling of the approximation with respect to 
the dimension of the SDE problem. 

§2.7| generalizes realizable discretizations in n-dimension to diffusion matrices 
that can be written as a sum of rank one matrices. 


• ^2.8 considers the special case of weakly diagonally dominant diffusion ma¬ 
trices, and proves that a second-order accurate, upwinded finite difference 
discretization is realizable on a gridded state space. 

Chapter applies realizable discretizations to SDE problems. 

• ^3.2| (resp. j ]3.5[ ) provides numerical evidence that the approximation is sta¬ 
ble (i.e. geometrically ergodic) and accurate in representing moments at finite 
times, mean first passage times, exit probabilities, and the stationary distri¬ 
bution of a one-dimensional process with additive (resp. multiplicative) noise. 

• §3.3| presents an asymptotic analysis of the mean-holding time of the approx¬ 
imation, which is used to validate some numerical results. 

• §3.6| shows that the approximation eliminates nonphysical moves in the Cox- 
Ingersoll-Ross model from mathematical finance, is accurate in representing 
the stationary distribution of this process, and can efficiently handle singular¬ 
ities in SDE problems by using adaptive mesh refinement. 

• §3.7| illustrates that the approximation accurately represents the spectrum of 
the generator in a variety of non-symmetric diffusion processes with additive 
noise including SDE problems with nonlinear drift fields, internal discontinu¬ 
ities, and singular drift fields. 


• (3.9 shows that the approximation on a gridded state space (with variable 
step sizes) is accurate in representing the stationary distribution of a two- 
dimensional process with multiplicative noise. 

• ( 3.10| applies the approximation to a Lotka-Volterra process from population 
dynamics, which shows that the approximation eliminates nonphysical moves, 
is suitable for long-time simulation, and can capture extinction. 

• §3.11| tests the approximation on a Brownian dynamics simulation of a colloidal 
cluster immersed in an unbounded solvent, and shows that the jump size can 
be set according to a characteristic length scale in an SDE problem, namely 
the Lennard-Jones radius. 

Chapter analyzes a realizable discretization on a gridded state space. 

• §4.1| precisely describes the theoretical context, and in particular, hypotheses 
on the SDE problem are given in Assumption |4.1| 
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§4.2| derives a sharp stochastic Lyapunov function for the SDE solution and 
the approximation, and uses this Lyapunov function to show that both the 
SDE solution and the approximation are geometrically ergodic with respect 
to an invariant probability measure. 

§4.3| solves a martingale problem associated to the approximation, justifies a 
Dynkins formula for a large class of observables of the approximation, and 
applies this solution to analyze basic properties of realizations of the approx¬ 
imation. 

j4.5 proves that the approximation is accurate with respect to finite-time 
conditional expectations and equilibrium expectations by using a continuous¬ 
time analog of the Talay-Tubaro expansion of the global error. 

Chapter analyzes a realizable discretization on a gridless state space, 
introduces a mollified generator that uses variable step sizes, 
shows that this generator defines a bounded operator with respect to 




Borel bounded functions and induces a Eeller semigroup. 

^5.3 proves that this generator is accurate with respect to the infinitesimal 
generator of the SDE. 

§5.4| shows that the associated approximation has a sharp stochastic Lyapunov 
function provided that the mollification is not too strong. 

^5.5 proves that the associated approximation is accurate with respect to 
finite-time conditional expectations and equilibrium expectations by using a 
continuous-time analog of the Talay-Tubaro expansion of the global error. 
Chapter analyzes tridiagonal, realizable discretizations on gridded state spaces. 

derives an explicit formula for the invariant density of the approximation, 
quantifies the accuracy of the stationary density of the approximation, 
derives an explicit formula for the exit probability (or committor function) 


§6T 


? 6.3 


of the approximation. 

• §6.4| derives an explicit formula for the mean first passage time of the approx¬ 
imation. 

Chapter concludes the paper. 
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CHAPTER 2 


Algorithms 


2.1. Realizability Condition 


Here we present a time-continuous approach to numerically solve an n-dimensional 
SDE with domain H C M". The starting point is the infinitesimal generator: 

(2.1.1) Lf{x) = Df{x)^^{x) +tTace{D‘^f{x)a{x)(T{x)'^) 
associated to the stochastic process Y that satisfies the SDE: 

(2.1.2) dY = n{Y)dt + V2a(Y)dW , y(0) S H 

where W is an n-dimensional Brownian motion and /x : H —)■ R” and cr : H —>■ 
are the drift and noise fields, respectively. (The test function / : K" —)■ R in (2.1.1) 


is assumed to be twice differentiable.) Let L* be the (formal) adjoint operator to 
L. Given an initial density p(0,a;), the Fokker-Planck equation: 


(2.1.3) 


dp 


describes the evolution of the probability density of Y (t). The time evolution of a 
conditional expectation of an observable is described by the Kolmogorov equation: 


(2.1.4) 


du 


with initial condition u{0,x). A weak solution to the SDE (2.1.2) can be obtained 
by solving (2.1.3) or ( |2.1.4 ). In this paper we consider solvers for (2.1.2) that exploit 
this connection. 

The essential idea is to approximate the infinitesimal generator of the SDE in 


(2.1.1) by a discrete-space generator of the form: 

K 


(2.1.5) 


Qf{x) = '^qix,yi{x)) (^fiViix)) - f{x)^ 


i=l 


where we have introduced a reaction rate function g : D x D —>■ R, and K reaction 
channels: yi{x) S Q for i = 1,..., K and for every x G Q. (This terminology comes 
from chemical kinetics 1311 .) Let h he a spatial step size parameter and p be a 
positive parameter, which sets the order of accuracy of the method. We require 
that the generator Q satisfies 

(Ql) pth-order accuracy: 

Qf{x) = Lf{x) + 0{hP) V X e H 

(Q2) realizability: 

q{x,y)>{) Mx.yeVi x^y 


9 
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Condition (Ql) is a basic requirement for any spatial discretization of a PDE. If 
(Ql) holds, and assuming numerically stability, then the approximation is pth-order 


accurate on finite-time intervals (see e.g. Theorem 4.5.1 for a precise statement). 


Condition (Q2) is not a standard requirement, and is related to the structure of 
the Kolmogorov equation. In words, (Q2) states that the weights in the finite 
differences appearing in (2.1.51 must be non-negative. We call an operator that 
satisfies this condition realizable. Without this requirement the approximation 
would be a standard numerical PDE solver and limited to low-dimensional systems 
by the curse of dimensionality. In contrast, when (Q2) holds the generator Q 
is realizable, in the sense that it induces a Markov jump process X which can 
be simulated exactly in time using the Stochastic Simulation Algorithm (SSA) 
[371 m [mi ig. This simulation does not require explicitly gridding up any 
part of state space, since its inputs - including rate functions and reaction channels 
- can be computed online. We stress that the channels in (2.1.5) can be picked 


so that the spatial displacement {yi — x) is bounded for i = 1,..., K, where as a 
reminder K is the ^ of channels. 

To emphasize that simulating X is simple, let us briefly recall how SSA works. 


Algorithm 2.1 (Stochastic Simulation Algorithm [37] 1. Given the current 
time t, the current state of the process X(t) = x, and the K jump rates evaluated 
at this state {q{x, yi{x))}^i, the algorithm outputs the state of the system A(f-|-i5t) 
at time t + St in two sub-steps. 


(Step 


(Step 


1 ): obtains St by generating an exponentially distributed random variable 
with parameter 


2 ) 


K 

Mx) = '^q{x,yi{x)) 

1=1 

: updates the state of the system by assuming that the process moves to 
state yi{x) with probability: 


P(A(t -h St) = yi{x) I X{t) = x) 


q{x,yi{x)) 

Xix) 


Algorithm suggests the following compact way to write the action of the 
generator Q in (2.1.5): 


(2.1.6) Qf{x) = Xix) E{f{x + ^x) - fix)) 

where the expectation is over the random n-vector whose state space is the set 

of all K reaction channels from the state x and whose probability distribution is 


]P(Cx = y^ix)) 


qjx^yijx)) 

Xix) 


l<i<K . 


To summarize, this paper proposes to (weakly) approximate the solution of a 
general SDE using a Markov process that can be simulated using Algorithm |2.1[ 
whose updating procedure entails generating a random time increment in (Step 1) 
and firing a reaction in (Step 2). By restricting the firings in (Step 2) to reaction 
channels that are a bounded distance away from the current state a;, the approxi¬ 
mation is numerically stable, a point we elaborate on in §4.2| Note that the mean 
time increment (or mean holding time) in (Step 1) is A(a:)“^. Thus, the spatial 
discretization of L determines the average amount of time the process spends at 
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any state. Morever, note that the process satisfies: X{s) = x for s G [t,t + 6t). 
The latter property implies that sample paths of the process are right-continuous in 
time with left limits (or cadlag). Also note that the process only changes its state 
by jumps, and hence, is a pure jump process. Since the time and state updates 
only depend on the current state X{to), and the holding time in (Step 1) is expo¬ 
nentially distributed, this approximation has the Markov property and is adapted 
to the filtration of the driving noise. Thus, the approximation is a Markov jump 
process. As we show in Chapter the Markov and non-anticipating property of 
the approximation greatly simplify the analysis. 

Next we discuss the state space of the approximation, and then derive schemes 
that meet the requirements (Ql) and (Q2). 

2.2. Gridded vs Gridless State Spaces 

We distinguish between two types of discretization based on the structure of 
the graph that represents the state space of the approximation. We say that a 
graph is connected if there is a finite path joining any pair of its vertices. 

Definition 2.2.1. A gridded (resp. gridless) state space consists of points 
which form a connected (resp. disconnected) graph. 

We stress that both state spaces permit capping the spatial step size of the approx¬ 
imation. Moreover, the SSA on either state space only requires local evaluation of 
the SDE coefficients. Gridded state spaces are appealing theoretically because a 
generator on such a space can often be represented as a matrix, and in low dimen¬ 
sional problems, the spectral properties of these matrices (e.g. low-lying eigenvalues, 
leading left eigenfunction, or a spectral decomposition) can be computed using a 
standard sparse matrix eigensolver. When the state space is unbounded this matrix 
is infinite. However, this matrix can still be approximated by a finite-dimensional 
matrix, and the validity of this approximation can be tested by making sure that 
numerically computed quantities like the spectrum do not depend on the trunca¬ 
tion order. We will return to this point when we consider concrete examples in 
Chapter]^ In order to be realizable, the matrix Q associated to the linear operator 
(5 on a gridded state space must satisfy the Q-matrix property. 

(2.2.1) Qij > 0 for all i ^ j and ^ = 0 . 

3 

Indeed, the Q-matiYx. property is satisfied by the matrix Q if and only if its associ¬ 
ated linear operator Q is realizable. Traditional finite difference and finite volume 
discretizations have been successfully used to design Q-matTYx. discretizations when 
the generator of the SDE is self-adjoint and the noise entering the SDE is additive, 
see [an film Horn [Toil [ 85 ]. However, Q-matrix discretizations on a gridded state 
space become trickier to design when the noise is multiplicative even in low dimen¬ 
sion. This issue motivates designing spatial discretizations on gridless state spaces. 
Let us elaborate on this point by considering realizable discretizations in one, two, 
and then many dimensions. With a slight abuse of notation, we may use the same 
symbol to denote an operator Q and its associated matrix. The meaning should be 
clear from the context. The most general realizable discretization is given in §2.7| 
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2.3. Realizable Discretizations in ID 


SDEs in one dimension possess a symmetry property, which we briefly recall 
since we expect the proposed schemes to accurately represent this structure. To 
prove this property, it is sufficient to assume that the noise coefficient is differ¬ 
entiable and satisfies ij{x) > 0 for all a; G 17, where 17 is the state space of the 
approximation. Under these assumptions the generator (2.1.1) can be written in 
conservative form as: 

(2.3.1) {Lf){x) = -^{M{x)v{x)f{x))' 


v{x) 

where M{x) = (a^x))"^ and i'{x) is given by: 


'(x) = 


M{x) 


exp 


M(s) 


ds 


for some a G 17. Similarly, the SDE ( 2.1.2[ ) can be written in self-adjoint form: 

(2.3.2) dY = -M{Y)U'{Y)dt + M’{Y)dt + V2a{Y)dW 

where U{x) is the free energy of v{x) defined as: 

U{x) = -\og{iy{x)) . 

Let (•, ■)^ denote the i^-weighted inner product defined as: 


if, 9 ) 1 ^= / f{x)g{x)v{x)dx 
Jn 

where the functions f{x) and g{x) are square integrable with respect to the measure 
v{x)dx. The generator L for scalar, elliptic diffusions is self-adjoint with respect to 
this inner product i.e. 


(2.3.3) {Lf,g), = {f,Lg), 

for smooth and compactly supported functions f{x) and g{x). Note that by taking 
g{x) = 1 for every x G 17, then {Lf, g)v = 0 for all smooth and compactly supported 
f(x), which implies that 

L*i'(x) = 0 . 

Thus, i'{x) is an invariant density of the SDE. Furthermore, if the function v{x) is 
integrable over 17, i.e. Z = iy(x)dx < 00 , then v{x)/Z is an invariant probability 
density (or stationary density) of the SDE. 


We stress that the generator of SDEs in more than one dimension often do not 
satisfy the self-adjoint property (2.3.3). Thus, to avoid losing generality, we will 
not assume the generator of the SDE is self-adjoint. In this context we present 
finite difference, finite volume, and probabilistic discretizations of L. 


2.3.1. Finite Difference. Let {xi] be a collection of grid points over 17 C IR 
as shown in Figure |2.1[ The distance between neighboring points is allowed to be 
variable. Thus, it helps to define forward, backward and average spatial step sizes: 

Sx~^ —)— Sx 

(2.3.4) 6 xf = Xi+i - x^ , 6 x~ = Xi - Xi_i , Sx^ = — - , resp. 

at every grid point Xi. The simplest way to construct a realizable discretization of 
L is to use an upwinded and a central finite difference approximation of the first 
and second-order derivatives in L, respectively. To derive this discretization, it is 
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helpful to write the generator L in (2.3.1) in terms of the functions M{x) and U{x) 
as follows: 


(2.3.5) 


{Lf){x) = -M{x)U'{x)f{x) + {M{x)f{x)y . 


To obtain a realizable discretization, the first term of ( |2.3.5 ) is discretized using a 
variable step size, first-order upwind method: 


5 xt 


6 x, 


Here we employ the shorthand notation: fi = f{xi) for all grid points Xi (and 
likewise for other functions), and 

a V 6 = max(a, b) and a Ab = min(a, b) 

for any real numbers a and b. By taking the finite difference approximation on the 
upstream side of the flow, i.e., going backward or forward in the finite difference ap¬ 
proximation of /' according to the sign of the coefficient —M^Uy this discretization 
is able to simultaneously satisfy first-order accuracy (Ql) and realizability (Q2). 
A variable step size, central scheme can be used to approximate the second-order 


term in (2.3.5): 




25xi 


(M,+i + - (M,_i + M,) 


Sx, 


It is straightforward to verify that this part of the discretization obeys (Ql) with 
p = 2 and since M{x) >0 satisfies (Q2) as well. The resulting discrete generator 
is a Q-matrix with off-diagonal entries 


(2.3.6) 


{Quk^-l = ^ A 0) + 


Note that to specify a Q-matrix, it suffices to specify its off-diagonal entries and 


then use (2.2.1) to determine its diagonal entries. Also, note that this scheme 


does not require computing the derivative of M{x), which appears in the drift 


coefhcient of (2.3.2). The overall accuracy of the scheme is first-order in space. 


The approximation also preserves an invariant density, see Proposition |6.1.1 for a 
precise statement and proof. Remarkably, the asymptotic mean holding time of 
this generator is exact when, roughly speaking, the drift is large compared to the 
noise. We will prove and numerically validate this statement in §3.3| 


2.3.2. Finite Volume. Next we present a finite volume discretization of L. 
This derivation illustrates the approach taken in |85] to construct realizable dis¬ 
cretizations for self-adjoint diffusions. For this purpose discretize the region H 
using the partition shown in Figure |2.2| This partition consists of a set of cells 


{Hi} with centers {xi\ spaced according to (2.3.4). The cells are defined as = 
ki-i/ 2 )Xi+i/ 2 \ where the midpoint between cell centers Xij^i /2 = {xi + Xi+i)/2 
gives the location of each cell edge. From (2.3.4) note that the length of each cell 
= 6 x,,. 


IS X. 


i+l/2 ~ Xi-1/2 
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fi-i 


h 



Xi-l 


Xi 


Xi-\-\ 


Figure 2.1. Finite Difference Grid. Consider non-uniformly 
spaced grid points {xj} and corresponding function values {fj = /(%)}■ 
This diagram depicts a three-point stencil centered at grid point Xi 
with corresponding function values. This stencil is used to approxi¬ 
mate the generator L using a realizable, finite difference approximation. 


See (2.3.61 for the entries of the resulting Q-matrix. 


Integrating the Fokker-Planck equation (2.1.3) over cell fli yields: 

d p{t,xy 




dx v(x) 


dx 


= ( u{x)M{x) 


d p{t,x) 
dx v{x) 


^i+l/2 


- 1/2 


— Pi+1/2 — Pi-l/2 

Here we have introduced ^.nd which are the probability fluxes at 

the right and left cell edge, respectively. Approximate the solution p{t, x) to the 
Fokker-Planck equation by its cell average: 

Pi{t) « 7 ^ / p{t,x)dx . 

Jili 

In terms of which, approximate the probability flux using 


Vi+l/2 ~ Mi+1/2 


(2.3.7) 


F, 


i-l-l/2 



Mi / Pi+i _ jH 

Sxf V^i-I-I Ui 
d p{t,x) 


.dx v{x) 

These approximations lead to a semi-discrete Fokker-Planck equation of the form: 

where Q* is the transpose of the matrix whose off-diagonal entries are 

1 Mi+i + Mi Ui+i -|- i^i 


(2.3.8) 



Sxi Sxf 2 2 vi 

1 Mi-i + Mi Vi-i + Vi 


Sxi 5 x, 


2vi 


This generator is a second-order accurate approximation of L and it preserves the 
invariant density v of the SDE in the following way. Let v be the invariant density 
V weighted by the length of each cell 

Vi = Sxi V. . 
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fi-i fi f' 


i+l 


Xi-1 


Xi 


Xi-\-\ 


Figure 2.2. Finite Volume Cells. Consider non-uniformly spaced 
points {xj} and corresponding function valnes {fj = f{xj)}. This dia¬ 
gram depicts a set of cells with cell centers given by the points {xj} and 
with cell edges located in between them at {xj +Xk)/‘2 for k = + 

This notation is used to derive a realizable, finite volume discretization 


of the generator L. See (2.3.81 for the resulting Q-matrix. 


It follows from (2.3.8) that 


^i+1 


since from (2.3.4) Sxf = This identity implies that u is an invariant density 


of see Proposition 


6 . 1.1 


for a detailed proof. Other approximations to the 


probability flux also have this property e.g. the midpoint approximation of Ui_|_i /2 
in (2.3.7) can be replaced by: 

Ui+i + Ui 


fo-ri/2 ~ fo A Uj+i or ^^+1/2 « exp 

(As a caveat, note that the former choice reduces the order of accuracy of the 
discretization because it uses the minimum function.) In fact, the evaluation of 
U can be bypassed by approximating the forward and backward potential energy 
differences appearing in (2.3.8) by the potential force to obtain the following Q- 
matrix approximation: 

5xt' 


(2.3.9) 



_ 1 -b Mi 

ic)i,i+l y r -I- 7, 

OXi oxj 2 

1 Mi i + Mi 


Sxi Sxi 


exp 


-U' 


St 

u' 

* 2 


In terms of function evaluations, (2.3.9) is a bit cheaper than (2.3.8), and this 
difference in cost increases with dimension. Under natural conditions on M(x) 
and U'{x), we will prove that (2.3.9) possesses a second-order accurate stationary 
density, see Proposition |6. 2.1 


For a more general theoretical and numerical treatment of realizable finite- 
volume methods for self-adjoint diffusions with additive noise we refer to 


2.3.3. Time-Continuous Milestoning. Finally, we consider a ‘probabilis¬ 
tic’ discretization of the generator L. The idea behind this Q-matrix approxima¬ 
tion is to determine the jump rates according to the exit probability and the mean 
holding time of the SDE solution at each grid point. 

To this end, let rjf = inf{t > 0 | F(t) = a, F(0) = x} be the first time the 
SDE solution hits the state ‘a’ given that it starts at the state ‘x’. Define the exit 
probability (or the committor function) to be q{x) = P(t^ < t“), which gives the 
probability that the process first hits ‘b’ before it hits ‘a’ starting from any state 
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‘x’ that is between them: a < x <b. This function satishes the following boundary 
value problem: 


(2.3.10) 


Lq{x) = 0 , q{a) = 0 , q{b) = 1 . 


Let Zi = q{xi) for every grid point Xi. Solving (2.3.10) with a = Xi-i and b = Xi+i 
yields the following exit probabilities conditional on the process starting at the grid 
point Xi 


(2.3.11) 




- Zj-I 
Zi+l — Zi — \ 
Zj+l - Zj 
Zi+1 ~ Zi_i 


We use these expressions below to define the transition rates of the approximation 
from Xi to its nearest neighbors. 

To obtain the mean holding time of the approximation, recall that the mean first 
passage time (MFPT) to (xi-i, Xi+i)^^ given that V(t) = x is the unique solution 
to the following boundary value problem: 


(2.3.12) Lu{x) = —1 , Xi-i < X < Xi^i , u(xi-i) = u(xi+i) = 0 . 


In particular, the MFPT to (a;i_i,x^+i)^ given that the process starts at the grid 
point Xi is given by 

A J = u(xi) . 

Since the process initiated at Xi spends on average u(xi) time units before it 
first hits a nearest neighbor of Xi, the exact mean holding time in state Xi is given 
by: 


(2.3.13) 




Similarly, the solution (2.3.11) specifies the transition probabilities: 
in.). 

(2.3.14) 


- Zi-i 


(Qe)i,i+1 “f -^i+l Zi—i 

-^ 2+1 Zi 


“f ■^z+1 1 


Incidentally, these are the transition probabilities used in discrete-time optimal 
milestoning [nni, where the milestones correspond to the grid points depicted in 


Figure 2.3 Combining (2.3.13) and (2.3.14) yields: 


(2.3.15) 




- 1-1 


Zj - Zj-i 1 

Zi+I - Zi-I u{Xi) 




Zj+i - Zj 1 
Zi+I - Z^-l u{Xi) 


which shows that a Q-matrix can be constructed in one-dimension based on the 
exact mean holding time and exit probability. We emphasize that this discretiza¬ 
tion is still an approximation, and in particular, does not exactly reproduce other 
statistics of the SDE solution as pointed out in ms]. We will return to this point 
in §6.4| where we analyze some properties of this discretization. In j |3.3| we will use 
this discretization to validate an asymptotic analysis of the mean holding time of 
finite-difference approximations in one dimension. 
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fi—l fi fi+1 


^i — 1 Xi Xi-\-\ 


Figure 2.3. Milestones. Consider a non-uniform (bounded or un¬ 
bounded) grid, a set of milestones {xj} on this grid, and corresponding 
function values {fj = f{xj)} as shown in this diagram. Using this 
notation a realizable discretization which retains the exact transition 
probabilities and mean holding times of the SDE is given in (2.3.151. 


To recap, the construction of Qe requires the exit probability (or committor 
function) to compute the transition probabilities, and the MFPT to compute the 
mean holding time. Obtaining these quantities entails solving the boundary value 
problems (2.3.10) and (2.3.12). Thus, it is nontrivial to generalize Qe to many 
dimensions and we stress our main purpose in constructing Qe is to benchmark 
other (more practical) approximations in one dimension. 


2.4. Realizable Discretizations in 2D 

Here we derive realizable hnite difference schemes for planar SDEs. The main 
issue in this part is the discretization of the mixed partial derivatives appearing in 


(2.1.1), subject to the constraint that the realizability condition (Q2) holds. To be 
specific, consider an SDE whose domain H is in and write its drift and noise 
fields as: 


Kx,y) = 




and a[x)a{x)^ = M{x,y) = 


M'^'^(x,y) 

M^^{x,y) 


M^‘^{x,y) 

M‘^\x,y) 


Here M is the diffusion field. Let S = {{xi^yj)} be a collection of grid points on 
a rectangular grid in two-space as illustrated in Figure 2.4 Given / : ; 

adopt the notation: 


we 


shorthand 

meaning 

fi,j ~ / (Xi, yj ) 

grid point function evaluation 

fi±i,j — fixi±i,yj) 
fi,j±i — f{xi,yj±i) 

nearest neighbor grid 
point function evaluation 


next to nearest neighbor 
grid point function evaluation 


The distance between neighboring vertical or horizontal grid points is permitted to 
be variable. Hence, it helps to define: 

5xf = Xi+i - Xi , 

ht = Vj+i - Vi > 
at all grid points {xi,yj) G S. 


Sx^ = Xi — Xi-i , Sxi = 


ox- +Ox, 


= yj - Vj-i . = 


Sy^ + ^y] 
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Figure 2.4. Grid in two dimensions. This diagram illustrates 
the reaction channels from the point (xi,yj) on a gridded space in two 
dimensions. These channels are used to obtain a realizable discretization 
of the generator L of the SDE. Note that the channels in the diagonal 
and antidiagonal directions are needed to approximate mixed partial 
derivatives that may appear in the continuous generator L. The result¬ 
ing realizable discretizations can be found in (2.4.61 and (2.4.71. 


We first consider a generalization of the first-order scheme (2.3.6). The drift 
term in L is approximated using an upwinded finite difference: 

(2.4.1) /i {xi,yj) ^^{xi,yj) ^ fi,j) Jx” /*—ij) 


and similarly, 


(mL-VO) 


(2.4.2) fj. ixi,yj) — {xi,yj) - 


+ - fi,j) + 


a 0) 

^y7 




We recruit all of the channels that appear in Figure |2.4| in order to discretize the 
second-order partial derivatives in L. Specifically, 

(2.4.3) 

trace{M(xi,yj)D'^f{xi,yj)) « - fij) + (/ij±i - fi,j) 

■' -^^ ^^-^ 

horizontal differences vertical differences 

+ C±(/,±i,j±i - /,,,) + (/ttij^i - /.,,) 


diagonal differences 


antidiagonal differences 


where we introduced the following labels for the rate functions: A^, and 

D^. To be clear there are eight rate functions and eight finite differences in this 
approximation. We determine the eight rate functions by expanding the horizon¬ 
tal/vertical finite differences in this approximation to obtain 


fi±l,j fi,j 
fi,j±l ~ fi,j 




Sx, 






A ± 


1 W 

2 \(9a;2 

1 W 

2 \dy'^ 


±\2 


) 






(2.4.4) 
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Similarly, we expand the diagonal/antidiagonal differences to obtain 


(2.4.5) 




k 


±1JT1 


-k. 


df 


dx 


± 5xf±\k- 


df 


dy 


1 

2 \dx'^ J 

(—) 

\dxdyj^ 


iSxff + 




SxfSy^ 




± 


dl 

dx 


SxfT 




dl 

dy 


2 ^9x2 J 

-( — ) 
\dxdyj^ 






SxfSy"^ 


Sy\ 


1 

2 W), 


,±\2 


(Sy,) 


Syf 




1 

2 WJ, 


{Syfl 


Note that the first order partial derivatives appearing in these expansions appear 
in plus and minus pairs, and ultimately cancel out. This cancellation is necessary 
to obtain an approximation to the second-order pure and mixed partial derivatives 
in (2.4.3). Using these expansions it is straightforward to show that the following 
choice of rate functions 

i± _ is/ri2 A n\ /'.sr..v.±A.i.T'i-i 


= MlkSxfSxl-^ - {Mk V 0) {SxfSyfr^ 
= M^^lSyfSy,)-^ - ^ q ) (SxfSyf)-^ - 

±1-1 


(Mi2 AO) (SxfSyf) 

rl2 


AO) (SxfSyf) 


Jl-i 
± 1-1 


V 0) (SxfSyf) 

D^ = -iMl^^A0) iSxfdyJ) 


imply that the approximation in (2.4.3) holds. Note here that we have used an 
upwinded finite difference to approximate the mixed partial derivatives of L. Com¬ 
bining the above approximations to the drift and diffusion terms yields: 


QufiJ — 


(2.4.6) 


mu 


+ - 


kj V 0 _ 

Sxf SxiSxt 

k.-i A 0 


M/2v0 M12 a0\^ 

' ^]ik+i,j-kj) 


Sxf Syl Sxf Sy^ 


Sx, 


+ 


mu 

^,3 




k,j V 0 

Syf 


SxiSx^ 
M22 

LI 

SyjSyf 


Ml2 V 0 
6 x-Sy- 

MU A0\ 


MU AO 
5x~Syf 


fi,j) 


MU VO 
SxlSyl 


+ 


fej dy~'' 


ikj+i - kk 


+ - 


k,j A 0 


-h 


\ 

Ml2 V 0 
Sx+Sy] 


+ 


M22 


MU VO 


SyjSy. Sx^ 5y. SxfSy, 


M12a0\, 

~ ) (/ij-l ~ fij) 


-ik 


MU VO 

+ lj + l ~ fij) + r - r - (/i-lj-1 ~ fij) 
dx^ dy. 


M/JaO^^ ^1 M/jAO^^ ^ ^ 

5xJdy, 5Xi 6yJ 
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This generator is a two-dimensional generalization of (2.3.51. However, unlike in 
ID, to approximate mixed partial derivatives finite differences in the diagonal and 
antidiagonal channel directions are used. To be accurate, the spurious pure deriva¬ 
tives produced by these differences, however, must be subtracted off the vertical 
and horizontal finite differences. This subtraction affects the realizability condition 
(Q2), and in particular, the weights in the horizontal and vertical finite differences 
may become negative. Note that this issue becomes a bit moot if the diffusion 
matrix M{x,y) is uniformly diagonally dominant, since this property implies that 

Ml] A m2] > I M]] I for all {x,,yj) G S 

and hence, Qu satisfies the realizability condition (Q2) on a gridded state space, 
with evenly spaced points in the x and y directions. 

Alternatively, if we approximate the first and second-order ‘pure’ partial deriva¬ 
tives in L using a weighted central difference, and use an upwinded finite difference 
to handle the mixed partials, we obtain the following generator: 

(2.4.7) 


Qcfi,j = exp ( +^yljSxf 


+ exp 


-I- exp 


-I- exp 


-I- exp 


-I- exp 


— exp 


exp 






l-2r- 


Ml] 

SxiSxf 

( Ml] 

ySxiSx^ 

( Mf] 

/ Ml] 

\ ^yj^yi 


Ml] VO 
Sxf6y] 


Ml] V 0 


Ml] A 0 


5x]5y. 


- fi,j) 


Sx^ Sy. Sx^ Sy 


Ml] A 0 


SxfSy] Sx^ Sy 


+-yh5xt -b -yljSyf 


Ml] V 0 

’ yt 

Ml] V 0 
6x-SyJ 
Ml] V 0 


M]] AO 


fij) 


ifi,j + l fi,j) 


Ml] A 0 


SxfSy 


- (/ij-l fi,j) 


SxfSy. 


(/, 


i-l-l J + l 


- kj) 


1 ~1 r - 1 ~2 


Sy', 


-nylMt - T^y-liSy'i 


-l^ljS^, +li^hSyj 


Mk VO, 

Sx, Syj 
Mil AO, 

^ + r - {fi+lj-l ~ fij) 
Sxl Sy^ 

Mil AO, ^ 

—1 „ + - fik ■ 


Sxi SyJ 

Here we have defined y = {k, kk the solution to the linear system 

M{x,y)fl{x,y) = kx,y) ■ 

Like Qu, if the diffusion matrix is dia gona lly d omin ant, then this generator is 
realizable on an evenly spaced grid. In (3.9 and { 3.10 we test Qc on planar SDE 
problems. In these problems the diffusion matrices are not diagonally dominant, 
and we must use variable step sizes in order to obtain a realizable discretization. 
In (2.8 we generalize the upwinded finite difference approximation Qc to multiple 
dimensions and prove that if the diffusion matrix is uniformly diagonally dominant, 
then this generalized discretization is realizable. 
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To avoid assuming the diffusion matrix is diagonally dominant, we next consider 
approximations whose state space is gridless. 

2.5. Realizable Discretizations in nD 

We now present a realization discretization on a gridless state space contained 
in C K". To introduce this approximation, set the number of channel directions 
equal to the dimension of the space and let i be an index over these reaction channel 
directions 1 <i <n. Let hf{x) be forward/backward spatial step sizes at state x, 
let hi{x) be the average step size at state x defined as 

h,{x) = {h+{x) + h-{x))/2 , 

let {ai{x)}2^i be the columns of the noise matrix cr(x), and let /l(x) be a trans¬ 
formed drift field defined as: 

(2.5.1) M{x)fl{x) = fi{x) 

where M{x) = a(x)(T(x)^ is the nxn diffusion matrix. (If M{x) is positive definite, 
then this linear system has a unique solution /i(a:) for every x S K". This condition 
holds in the interior of fl for all of the SDE problems we consider in this paper.) 
With this notation in hand, here is a second order accurate, realizable discretization 
in n-dimensions. 


(2.5.2) 



The generator Qc is a weighted central finite difference scheme in the directions 
of the columns of the noise coefficient matrix evaluated at the state x. When 
the noise coefficient matrix is additive (state-independent) and isotropic (direction- 
independent), then the state space of this generator is gridded and the generator 
can be represented as matrix. In general, we stress that the state space of this 
approximation is gridless, and thus the generator cannot be represented as a ma¬ 
trix. Nevertheless it can be realized by using the SSA algorithm. A realization of 
this process is a continuous-time random walk with a user-defined jump size, jump 
directions given by the columns of the noise coefficient matrix, and a mean holding 
time that is determined by the discretized generator. This generator extends real¬ 
izable spatial discretizations appearing in previous work IMlITim ITO IT041I85] 
to the SDE (2.1.2), which generically has a non-symmetric infinitesimal generator 


L and multiplicative noise. 

The next proposition states that this generator satisfies (Ql) and (Q2). 
Proposition 2.5.1. Assuming that: 

~ \ all x G 

then the operator Qc in (2.5.2[) satisfies (Q2) and (Ql) with p = 2. 
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Proof. This proof is somewhat informal. A formal proof requires hypotheses 
on the SDE and the test function f{x). For clarity of presentation, we defer these 
technical issues to Chapter]^ Here we concentrate on showing in simple terms that 
the approximation meets the requirements (Ql) and (Q2). 

Let D^f{x){ui, ..., Uk) denote the /cth-derivative of f{x) applied to the vectors 
ui,...,Uk for any fc > 1. Note that the rate functions appearing in Qc are all 
positive, and thus, the realizability condition (Q2) is satisfied. To check second- 
order accuracy or (Ql) with p = 2, Taylor expand the forward and backward finite 
differences appearing in the function Qcf(x) to obtain: 


Qcf{x) = ^ j^^Df{x){a,{x)) ^exp 

-exp 

+ f(.x){a^{x),a^{x)) (^h+{x)exp (^fl{xf 

+h~{x)exp 

+ -^j:^D^f{x){cr^{x),a^{x),a,{x)) (^hf {xf exp 

-h-{xfexp + 0{h^) . 


In order to obtain a crude estimate of the remainder term in this expansion, we 
used the identity -|- 6^ = (a -I- b){a? + — ab) with a = h^l {x) and b = h~{x), 

and the relation hi{x) = {hf(x) + h~{x))/2. A similar Taylor expansion of the 
exponentials yields. 


Qcf{x) = ^ Df{xY'ai{x)il{xY'(Ti{x) + D^f{x){ai{x), ai{x)) + 0{h^) 
2 = 1 


= trace 


n 

(^Df{x)fi{xf -1- D^f{x)^ CTiCrf 


2 = 1 


+ 0(h2 


=M(x) / 


= trace {^Df (x)p,{x)^ + f{x)M{x)) + 0{h^) 

= Lf{x) + 0{h^) 


where we used the cyclic property of the trace operator, the basic linear algebra 
fact that the product trcr^ is the sum of the outer products of the columns of a, and 
elementary identities to estimate the remainder terms, like — b'^ = {a + b)(a — b) 
with a = hf{x) and b = h~{x). Thus, Qc is second-order accurate as claimed. □ 


The analysis in Chapter focusses on the stability, accuracy and complexity 
of the generator Q^, mainly because it is second-order accurate. However, let us 
briefly consider two simpler alternatives, which are first-order accurate. To avoid 
evaluating the exponential function appearing in the transition rates of Qc in (2.5.2), 
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an upwinded generator is useful 


Qufix) = '^2 ^ ^ if i^ + ix)criix)) - f{x)) 


h+ix) 

{pL^x)'^ai{x)) A 0 


(2.5.3) 


if{x-hi {x)a^{x)) - f{x)) 


K {x) 

,+( ^ if {x)(r^{x)) - fix)) 

hjix)hiix) 

if ix - K ix)a^ix)) - fix)) . 


h- ix)hiix) 

Let be the standard basis of K”. To avoid solving the linear system in 

(2.5.1|), an upwinded generator that uses the standard basis to discretize the drift 


term in L is useful 


Quufix) = '^2 ^ ^ if i^ + (2;)e*) - fix)) 


K (x) 

(/r(a;)^ei) A 0 


(2.5.4) 


K ix) 

1 


if ix - K ixy^) - fix)) 


hfix)hiix] 

1 


if ix + hf ix)aiix)) - fix)) 
if ix - K ix)aiix)) - fix)) . 


ix)hiix) 

Proposition 2.5.2. Assuming that: 

\hfix) — h~ix) \ V hfix) V h~ix) < h for all x € ft 
then the operators Qu and Quu satisfy (Q2) and (Ql) with p = 1. 


The proof of Proposition |2.5.2| is similar to the proof of Proposition |2.5.1| and 
therefore omitted. 

2.6. Scaling of Approximation with System Size 

Here we briefly discuss the complexity of the SSA induced by the generator in 
(2.5.21. We will return to this point in Chapter]^ where we quantify the complexity 
of sample paths as a function of the spatial step size parameter. 

First we compare the complexity of the SSA to a standard numerical PDF 
solver for the Kolmogorov equation. When the grid is bounded, even though the 
matrix associated to standard PDF discretizations (like finite difference, volume and 
element methods) is often sparse, the computational effort required to construct its 
entries in order to time integrate the resulting semi-discrete ODF (either exactly 
or approximately) scales exponentially with the number of dimensions. On an 
unbounded grid, this ODF is infinite-dimensional and a standard numerical PDF 
approach does not apply. 


In contrast, note from Algorithm 2.1 that an SSA simulation of the Markov 
process induced by the generator in (2.5.21 takes as input the reaction rates and 


channels, which can be peeled off of (2.5.2). More to the point, both the spatial and 


temporal updates in each SSA step use local information. (Put another way, in the 
case of a gridded state space, a step of the SSA method only requires computing a 
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row of a typically sparse matrix associated to the generator, not the entire matrix.) 
Thus, the cost of each SSA step in high dimension is much less severe than standard 
numerical PDE solvers and remains practical. 

Next consider the average number of SSA steps required to reach a finite time. 
Referring to Algorithm |2.1[ the random time increment ti — to gives the amount 
of time the process spends in state X(to) before jumping to a new state X{ti). In 
each step of the SSA method, only a single degree of freedom is updated because 
the process moves from the current position to a new position in the direction of 
only one of the n columns of the noise matrix. In other words, the method involves 
local moves, and in this way, it is related to splitting methods for SDEs and ODEs 
based on local moves, for example the per-particle splitting integrator proposed for 
Metropolis-corrected Langevin integrators in |14j or J-splitting schemes for Hamil¬ 
tonian systems [33] . However, the physical clock in splitting algorithms does not 
tick with each local move as it does in (Step I) of the SSA method. Instead the 
physical time in a splitting scheme does not pass until all of the local moves have 
happened. In contrast, the physical clock ticks with each local move of the SSA 
method. This difference implies that the mean holding time of the SSA method is 
inversely proportional to dimension. 

The easiest way to confirm this point is to consider a concrete example, like 


evolving n decoupled oscillators using the SSA induced by (2.5.2). In this example, 


the mean time step for a global move of the system - where each oscillator has 
taken a step - cannot depend on dimension because the oscillators are decoupled. 
For simplicity assume the SDE for each oscillator is dY = fj,{Y)dt + y^dW and the 
spatial step size in every reaction channel direction is uniform and equal to h, and 
let X = {xi,... ,Xn)'^- Then, from (2.5.2) the total jump rate is: 




i=l 


where Xi{x) is the jump rate of the *th oscillator: 

Xi{x) = p cosh(/r(a;i)^) . 
Expanding the total jump rate about h = 0 yields 


cy 1 ^ 

= IT + 7 H ■ 


i=l 


Similarly, expanding the mean holding time about h = 0 yields 






Thus, the total jump rate is roughly 0{{2n)/h?) and the mean holding time scales 
like 0{h?/{2n)) with dimension n, as expected. 

To summarize, we see that given the nature of the SSA method, accuracy 
necessitates a smaller mean time step for each local move and the overall complexity 
of the SSA method is no different from a splitting temporal integrator that uses 
local moves. 
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2.7. Generalization of Realizable Discretizations in nD 

Suppose the diffusion matrix can be written as a sum of rank-1 matrices: 

S 

(2.7.1) M{x) = y^Wi{x)'qi{x)rii{xY' 

i=l 


where rjiix) S K” is an n-vector and Wi{x) > 0 is a non-negative weight function for 
all x G n, and for 1 < i < s. Since the diffusion matrix M{x) is positive definite, 
it must be that s > n and that the set of vectors spans M" for every 

X G M". For example, since M{x) = a{x)a{x)’^, the diffusion matrix can always be 
put in the form of ( 2.7.1| ) with s = n, unit weights Wi{x) = 1, and rii{x) = ai{x) 
for 7 = 1,..., n, or equivalently. 


M{x) ='^Ui{x)ai{xY' . 

i=l 


The following generator generalizes (2.5.2) to diffusion matrices of the form (2.7.1). 


(2.7.2) 



This generator is a weighted, central finite difference approximation along the di¬ 
rections given by the vectors {r]i{x)}f^i. To be sure /i in (2.7.2) is the transformed 


drift held introduced earlier in ( 2.5.1| ). The next proposition states that this gen¬ 
eralization of Qc meets the requirements (Ql) and (Q2). 


Proposition 2.7.1. Suppose that the diffusion matrix satisfies (2.7.1). 
suming that: 


As- 


\hf{x) — /ij (a;)| V {hf{x)Y V {h^ (x))^ < hfi for all x gVI 


then the operator Qc in (2.7.2) satisfies (Q2) and (Ql) with p = 2. 


Proof. The operator Qc satishes the realizability condition (Q2) since wfix) > 
0 by hypothesis on the diffusion held M{x). To check accuracy, Taylor expand the 
forward/backward diherences appearing in Qcf{x) to obtain: 


Qcf{x) = ^ ^exp - exp 

+ ^0"^{ht exp + K exp 

+ - (h)-)^ exp + 0{h^) . 
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A similar Taylor expansion of the exponentials yields, 


Qcf{x) = ^ w,Df'^r]ijFr]i + WiD^f{r]i, 77 ,) + 


/ 


= trace 


+ D^f) 


( \ 


V 




0{h^ 


V =M(x) / 

= trace (^Df{x) + D^f{x) M{x)) + 0{h^) 

= Lf{x) + 0{h'^) 


where we used (2.7.1). Thus, Qc is second-order accurate as claimed. 


□ 


As an application of Proposition |2.7.1| we next show that realizable discretiza¬ 
tions on gridded state spaces are possible to construct if the diffusion matrix of the 
SDE is diagonally dominant. 


2.8. Weakly Diagonally Dominant Case 


Here we extend the upwinded finite-difference approximation Qc in (2.4.7) from 


two to many dimensional SDEs assuming that the diffusion matrix is weakly di¬ 
agonally dominant. We prove that this extension is realizable and second-order 
accurate under this assumption. Let us quickly recall what it means for a matrix 
to be weakly diagonally dominant. 


Definition 2.8.1. An n x n matrix A is weakly diagonally dominant if there 
exists an n-vector y = (?/i,..., with positive entries such that 


( 2 . 8 . 1 ) 


\ Aij\yj 


for alH € {1,..., n}. 

The following lemma is straightforward to derive from this definition. 


Lemma 2.8.1. Let A be a square matrix satisfying (2.8.1). Then there exists 


a positive diagonal matrix P such that the square matrix A = PAP is diagonally 
dominant. 

Let A(a:) = (Ai(a;),..., A„(a;))'^ be a field of spatial step sizes, which permit 
a different spatial step size for each coordinate axis. Let fl be the transformed drift 


field defined in (2.5.1). Here is a multi-dimensional generalization of Qc 


(2.8.2) 


Qcf{x) = \{x) E(/(a; -h fx) - f{x)) 
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where the expectation is over the random vector which has a discrete probability 
distribution given by: 


^ HxY 


f-SOI 

yi=l ^ V * 

+ ^ <^( 2 /z”)exp -a;)") 


2=1 






(2.8.3) 


+ E 


l<2<n 


j>i 


+ E 

Hy-/) 

l<2<n 


j>i 


1 

M 

S{ytn 

l<2<n 


j>i 


1 

M 

Siy-fl 


exp 


- x) 




l<2<r 

j>i 


Mij V 0 

AiA. 


M,j V 0 

A.A, 


A 0 

A,Aj 


Mij A 0 
A,A, 


where <5(1/) is a point mass concentrated at y, X{x) is a normalization constant, 
A(a:)“^ is the mean holding time, and we have introduced the following reaction 
channels: 


(2.8.4) 


yf = x± CiAi , 1 < z < n 

±,d 


= X ± e^A.i ± CjAj , l<i<7i,j>z 

jL-ij , 1 < z < n , j > 1 


Vij 

= a; ± CiAi =f e,-A^ 


The superscripts d and a in the reaction channels above refer to the diagonal and 
anti-diagonal directions in the i_/-plane, respectively. 

The following proposition states that Qc satisfies (Ql) and (Q2). 


Proposition 2.8.2. Consider the generator Qc given in (2.8.2). If for all 
X € M", the Cholesky factorization of the diffusion matrix M(x) has at most two 


nonzero entries in every column, then there exists a field A such that Qc in (2.8.2) 
is realizable and second-order accurate. 


Proof. This proof is an application of Proposition |2.7.1| To apply this propo¬ 
sition we show that if the Cholesky factorization of M (x) has at most two nonzero 


entries per column, then M{x) can be put in the form of (2.7.1). According to 


Theorem 9 of [8], the Cholesky factorization of a symmetric matrix with positive 
diagonals has at most two nonzero entries in every column if and only if it is weakly 
diagonally dominant. Since M(x) is symmetric positive definite by hypothesis, and 
since the Cholesky factorization of M(x) has at most two nonzero entries per col¬ 
umn also by hypothesis, M(x) is weakly diagonally dominant for all x € K". Hence, 
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by Lemma 2.8.1 there exists a positive diagonal matrix P(x) such that: 

(2.8.5) M{x) = P{x)M{x)P{x) 

is diagonally dominant for each x G K". Define a field of step sizes A by setting 
Ai = h/Pii. Hence, we have that: 

n 

M Maeicf + ^ MijideJ + e^-ef) 


2=1 


j>i 


ivt - ivt - ^ 


i=l j/i 

l<2<n 

l<2<n 

j>i 


where we used the forward direction channels defined in (2.8.4). Note that the 
coefficient of each rank one matrix in this last expression is non-negative because 
M{x) is diagonally dominant. Thus, M is in the form of (2.7.1) with s = and 
note that Qc is a special case of Qc in (2.7.2) with hf{x) = h. The desired result 
follows from Proposition |2. 7.1] □ 

We have the following corollary. 

Corollary 2.8.3. If the n x n diffusion matrix M{x) is diagonally dominant 
for all x € K", then Qc in (2.8.2) is realizable on a gridded state space in K". 

As far as we can tell, realizable, finite difference discretizations of this type 
seem to have been first developed by H. Kushner |77l [751 179]. 

Proof. In this case the diffusion matrix can be written as: 

n 




i=l 




T 'y ' 0) (Si T T 


l<i<n 

j>i 


- X! ^ 0) (ei - ej)(ei - 


l<i<n 

j>i 

and the coefficient of each rank one matrix in this last expression is non-negative 
because M{x) is diagonally dominant by hypothesis. Thus, M is in the form of 
( 2.7.1| ) with s = and note that Qc is a special case of Qc in (2.7.2) with hf{x) = 
h. Note also that the state space is gridded since the jumps are either to nearest 
neighbor grid points: x ± hei for all 1 < i < n, or to next to nearest neighbor grid 
points: x ± hei ± hej and x ± hci hcj for all 1 < * < n and j > i. □ 
















CHAPTER 3 


Examples &; Applications 


a cubic oscillator in ID with additive noise as described in ^3.2 


a log-normal process in ID with multiplicative noise as described in ^3.5 


the Cox-Ingersoll-Ross process in ID with multiplicative noise as described 
in 


3.1. Introduction 

In this chapter we apply realizable discretizations to 
( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 

( 9 ) 


a non-symmetric Ornstein-Uhlenbeck process as described in S 3.7.1 


the Maier-Stein SDE in 2D with additive noise as described in S 3.7.2 


simulation of a particle in a planar square-well potential well in 2D with 
additive noise as described in §3.7.3 


simulation of a particle in a worm-like chain potential well in 2D with 
additive noise as described in §3.7.4[ 


a log-normal process in 2D with multiplicative noise as described in { 3.9 


a Lotka-Volterra process in 2D with multiplicative noise as described in 
§3.10| and, 

(10) a Brownian dynamics simulation of the collapse of a colloidal cluster in 
39D with and without hydrodynamic interactions as described in §3.11| 
For the ID SDE problems, we numerically test two realizable spatial discretiza¬ 
tions with a gridded state space. The first uses an upwinded (resp. central) finite 
difference method to approximate the first (resp. second) order derivative in (2.1.1). 


(3.1.1) 


{Qu)i,i+1 — V 0) -I- 

OX, 


-{fi, A 0) 


M, 
6x. 
M, 

6xi 


The second generator uses a weighted central finite difference scheme to approxi¬ 


mate the derivatives in (2.1.1). 


(3.1.2) 


{Qc)i,i+1 — 




Mi exp 


^ 6xf \ 

M, 2 J 


. , I Mi 
— Mi exp — —— 
5xi Sx, \ Mi 2 


1 

5xi 5xf 

1 


These discretizations are slight modifications to the discretizations given in (2.3.6) 


and (2.3.9), respectively. For the planar SDE problems, we apply a 2D version of 


(3.1.2) given in (2.4.7). When the noise is additive, we use an evenly spaced grid. 


When the noise is multiplicative, we use adaptive mesh refinement as described in 
§3.4|in ID and §3.8|in 2D. For the 39D SDE problem, we apply the generator on a 


gridless state space given in (2.5.2). 
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3. EXAMPLES & APPLICATIONS 


3.2. Cubic Oscillator in ID with Additive Noise 

This SDE problem is a concrete example of exploding numerical trajectories. 


Consider (2.1.2) with = M, ii{x) = —x^, and cr{x) = 1 i.e. 
(3.2.1) dY =-Y^dt + V2dW , r(o)eK. 


The solution to (3.2.1) is geometrically ergodic with respect to a stationary distri¬ 


bution with density 
(3.2.2) p(a;) = exp(— a;'^/4) , where Z = f exp{—x^/A)dx. 

Jk 

However, explicit methods are transient in this example since the drift entering 
( 3.2.1| ) is only locally Lipschitz continuous. More concretely, let {Xn} denote a 
(discrete-time) Markov chain produced by forward Euler with time step size hf. 
Then, for any ht > 0 this Markov chain satisfies 


^xX'^t/ht\ 


oo as t —?> oo 


where denotes expectation conditional on Xq = x. (See Lemma 6.3 of |53j for a 
simple proof, and [60] for a generalization.) Metropolizing explicit integrators, like 
forward Euler, mitigates this divergence. Indeed, the (discrete-time) Markov chain 
produced by a Metropolis integrator {Xn} is ergodic with respect to the stationary 
distribution of the SDE, and hence. 




x‘^v{x)dx as t —)■ oo . 


However, the rate of convergence of this chain to equilibrium is not geometric, 
even though the SDE solution is geometrically ergodic The severity of 

this lack of geometric ergodicity was recently studied in [U, where it was shown 
that the deviation from geometric ergodicity is exponentially small in the time step 
size. This result, however, does not prevent the chain from “getting stuck” (rarely 


accepting proposal moves) in the tails of the stationary density in (3.2.2). 

In this context we test the generators Qu in ( |3.I.l ) and Qc in (3.1.2) on an in¬ 
finite, evenly spaced grid on M. Figure |3T] plots sample paths produced by the SSA 
induced by these generators. The initial condition is large and positive. For compar¬ 
ison, we plot a sample path (in light grey) of a Metropolis integrator with the same 
initial condition. At this initial condition, the Metropolis integrator rarely accepts 
proposal moves and gets stuck for the duration of the simulation. In contrast, the 
figure shows that the proposed approximations do not get stuck, which is consistent 
with Theorem |4.2.8| The inset in the figure illustrates that the time lag between 
jumps in the SSA method is small (resp. moderate) at the tails (resp. middle) of 
the stationary density. This numerical result manifests that the mean holding time 
adapts to the size of the drift. Note that the time lag for Qf. is smaller than for Qu¬ 
in the next section, this difference in time lags is theoretically accounted for by an 
asymptotic analysis of the mean holding time. Figure |3.2| provides evidence that 
the approximations are accurate in representing the stationary distribution, mean 
first passage time, and committor function. To produce this hgure, we use results 
described in Chapter In particular, in the scalar case, the numerical station¬ 
ary density, mean first passage time and committor satisfy a three-term recurrence 
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relation, which can be exactly solved as detailed in Sections |6.1[ |6.4| and |6.3[ re¬ 
spectively. We also checked that these solutions agree - up to statistical error - 
with empirical estimates obtained by an SSA simulation. 

Remark 3.2.1. It is known that the rate of convergence (and hence accuracy) 
of the solution of standard time integrators to first exit problems drops because they 
fail to account for the probability that the system exits in between each discrete time 
step [38l[3mi40) . Building this capability into time integrators may require solving 
a boundary value problem per step, which is prohibitive to do in many dimensions 
[96]. By keeping time continuous, and provided that the grid is adjusted to fit the 
boundaries of the region of the first exit problem, the proposed approximations take 
these probabilities into account. As a consequence, in the middle panel of Figure [3^ 
we observe that the order of accuracy of the approximations with respect to first 
exit times is identical to the accuracy of the generator, that is, there is no drop of 
accuracy as happens with time discretizations. To be precise, this figure illustrates 
the accuracy in these approximations with respect to the mean first passage time of 
the cubic oscillator to (0, 2)‘^ and shows that the generators Qu and Qc are 0{6x) 
and O(Sx^), respectively. 


3.3. Asymptotic Analysis of Mean Holding Time 


Consider the following instance of (2.1.2) 

(3.3.1) dY = fi(Y)dt + V2dW , y(0)eK. 

Assume the drift entering ( |3.3.1 ) is differentiable and satisfies a dissipativity con¬ 
dition, like 

sign(x)p-(a;) —)■ —oo as |a;| —>■ oo . 

Hence, for sufficiently large |a;| the SDE dynamics is dominated by the drift, and 
asymptotically, the time it takes for the SDE solution Y(t) to move a fixed distance 
in space can be derived from analyzing the ODE 


(3.3.2) 


Y = My), r(o) = x. 


Equation (3.3.2) implies that the time lapse between T(0) = Xi and Y(t*) = Xi-i 
where Xi > Xi-i ^ 0 satisfies 


(3.3.3) 


= 


dx 

._i \dix)\ ' 


By using integration by parts and the mean value theorem, observe that (3.3.3) can 
be written as: 


= 


Sx 


\Kx^)\ 

= t* - 


{x — Xi-i)fj,{x) fj,'{x)dx 




' 2^{(r 

for some ^ G [xi-i,Xi], and where we have introduced: t* = dx/\^i\. Let us compare 
t® to the mean holding times predicted by the generators: Qu in (3.1.1) and Qc in 
(3.1.2). For clarity, we assume that 6x = Sxf = Sx~. 
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Sample Paths 


10 ^ 

8 

6 

^ 4 



MALA 

-Qu-SSA 

—Qc-SSA 



Figure 3.1. Cubic Oscillator. The main figure plots sample paths 
produced by a Metropolis integrator (MALA) and the SSA induced by 
Qu in (2.3.61 and Qc in (2.3.91 over the time interval [0,100]. The tem¬ 


poral step size of MALA and the spatial step size of the SSA integrators 
are set at 5t — 0.015 and 5x = 0.25, respectively. The mean holding 
time for the SSA integrators operated at this spatial step size is ap¬ 
proximately (St) = 0.03. At the initial condition marked by the black 
dot at t = 0, the trne dynamics is dominated by the deterministic drift, 
which causes the true solution to drop toward the origin. Notice that 
MALA gets stuck at high energy, which illustrates that it is not geomet¬ 
rically ergodic. In contrast, SSA is geometrically ergodic (on an infinite 
grid), and hence, does not get stuck. The inset shows that the SSA time 
adapts to the size of the drift: when the drift is large (resp. small) the 
lag times are smaller (resp. larger). 


By hypothesis, = ^[xi) is less than zero if Xi is large and positive, and hence, 

oe I 

= {{Qu)i,i+1 + iQu)i,i-l) ^ = 


from (3.1.1 1 the mean holding time of Qu can be written as: 

Sx'^ 

2 -I- \fii\Sx 


This expression can be rewritten as 
(3.3.5) e = t* - f 


2 -b \ni\Sx 


From (3.3.51, we see that approaches t*, as the next Proposition states. 

Proposition 3.3.1. For any Sx > 0, the mean holding time of Qu satisfies: 

n I I 

-(J as X,; —>■ oo . 
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Likewise, if the second term in (3.3.4) decays faster than the first term, then the 
relative error between and t* also tends to zero, and thus, the estimate predicted 
by Qu for the mean holding time asymptotically agrees with the exact mean holding 
time. This convergence happens if, e.g., the leading order term in is of the 
form —a for p > 0 and a > 0. 

Repeating these steps for the mean holding time predicted by Qc yields: 


(3.3.6) 


Sx^ 

^ — ((l5c)z,i+i ^ sech(p2*^x) 


It follows from this expression that even though Qc is a second-order accurate 
approximation to L, it does not capture the right asymptotic mean holding time, 
as the next Proposition states. 


Proposition 3.3.2. For any Sx > 0, the mean holding time of Qc satisfies: 

1 as Iced —>■ oo . 


\F - t* 
t^ 


converges to zero too fast. This analysis 
which shows that the mean holding time for agrees 
with the mean holding time of the generator Qe that uses the exact mean holding 


Simply put, the mean hold ing time of Qc 
is confirmed in Figure 


3.3 


time. This generator was introduced in (2.3.3 see (2.3.15) for its definition. This 


numerical result agrees with the preceding asymptotic analysis. In contrast, note 
that the mean holding time of Qc is asymptotically an underestimate of the exact 


mean holding time, as predicted by (3.3.6). 


3.4. Adaptive Mesh Refinement in ID 

The subsequent scalar SDE problems have a domain Id = M+, and may have 
a singularity at the origin for certain parameter values. To efficiently resolve this 
singularity, we will use adaptive mesh refinement, and specifically, a variable step 
size infinite grid S = {xi} € K“'" that contains more points near the singularity at 
X = 0. We will construct this grid by mapping the grid S to logspace to obtain a 
transformed grid S' on M defined as 

S = = log(a:i) , y x^G S . 

Note that as —>■ —oo (resp. —>■ oo) we have that Xi ^ 0 (resp. Xi —>■ oo). For 
the transformed grid in logspace, we assume that the distance between neighboring 
grid points is fixed and given by <5^, i.e., 

for alH S Z . 


Since 


x^+i = exp(^j+i) = exp(i5^) exp(^i) = exp{S£,)xi 


it follows from the definitions introduced in (2.3.4) that: 


(3.4.1) Sxf = (exp(5^) — l)xi , 6x^ = (1 — exp(—i5^))xi , 6xi = sinh(5^)a:i . 
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3.5. Log-normal Process in ID with Multiplicative Noise 


Consider (2.1.21 with 


fl = , ^(cc) = —X log X + x , M{x) = {a{x))"^ = , 

and initial condition Y (0) £ fl. This process has a lognormal stationary distribution 
with probability density: 


(3.5.1) 


^ - log(2^) ) • 


In fact, log(y(t)) satisfies an Ornstein-Uhlenbeck equation with initial condition 
log(y(0)). Exponentiating the solution to this Ornstein-Uhlenbeck equation yields 

y (t) = exp log y (0)-I- J . 

It follows that: 

(3.5.2) E 2 ,(y(t)^) = exp (2e“‘log(a:)-I-2(1 — e“^*)) . 

The formulas (3.5.1) and (|3 .5.2[ ) are useful to numerically validate the approxima¬ 
tions. 


Figure 3.4 show numerical results of applying the generators in (3.1.1) and 


(3.1.2) to solve the boundary value problems associated to the mean first passage 


time and exit probability (or committor function) on the interval [1/2, 5] and using 
an evenly spaced mesh. Figure |3.5| illustrates that these generators are accurate 


with respect to the stationary density (3.5.1) using the adaptive mesh refinement 
described in §3.4[ Figure [T6| plots sample paths produced by SSA for these gen¬ 
erators and the SSA time as a function of the step index. Figure |3.7| illustrates 


the accuracy of the proposed approximations with respect to the mean-squared 
displacement Ea;(y(l)^) with an initial condition of y(0) = a; = 2. The relative 
error is plotted as both a function of the spatial step size (top panel) and the mean 


lag time (bottom panel). Finally, Figure 3.8 plots the average number of compu¬ 


tational steps (top panel) and the mean lag time (bottom panel) as a function of 
the spatial step size. The statistics are obtained by averaging this data over 100 
SSA realizations with initial condition X(0) = 1 and time interval of simulation of 
[ 0 , 10 ]. 


3.6. Cox-Ingersoll-Ross Process in ID with Multiplicative Noise 

Consider a scalar SDE problem with boundary 

(3.6.1) dY = P{a-Y)dt + aVYdW , y(0) G K+ . 

In quantitative finance, this example is a simple model for short-term interest rates 
[23] . It is known that discrete-time integrators may not satisfy the non-negativity 
condition: y(t) > 0 for all t > 0 [891 [5]. Corrections to these integrators that 
enforce this condition tend to be involved, specific to this example, or lower the 
accuracy of the approximation. This issue motivates applying realizable discretiza¬ 
tions to this problem, which can preserve non-negativity by construction. 
According to the Feller criteria, the boundary at zero is classihed as: 

• natural if 2/3a/cr^ > 1; 

• regular if 1 > 2/3a/(T^ > 0; and, 

• absorbing if j3a < 0. 
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When the boundary at zero is natural (resp. regular), the process is positive (resp. non¬ 
negative). When the origin is attainable, the boundary condition at the origin is 
treated as reflective. If /3 > 0 and a > 0, the SDE solution is ergodic with respect 
to a stationary distribution with density proportional to: 


(3.6.2) 


i'{x) = X 


i-i / 2/3 

exp [ 


We emphasize that v{x) is integrable over IR+ if and only if /3 > 0 and a > 0. 


Note that when the boundary at zero is regular, the stationary density in (3.6.2 1 is 
unbounded at the origin. 

The numerical tests assess how well the approximation captures the true long¬ 
time behavior of the SDE. To carry out this assessment, we estimate the difference 
between the numerical stationary densities and 


(3.6.3) 


y = Z- 




v{x)dx for all Xi G S 


where the constant Z is chosen so that — 1- Figures 3.9 and 3.10 illustrate 


the results of using adaptive mesh refinement. The top panels of Figures [3f^ and 


3.10 plot the free energy of the stationary density when the boundary at zero 
is natural and regular, respectively. The benchmark solution is the cell-averaged 


stationary density (3.6.2) on the adaptive mesh shown in the inset. Note that when 


the boundary at zero is natural, the minimum of the free energy of the stationary 
density is centered away from zero. In contrast, when the boundary is regular, the 
minimum of the free energy is attained at the grid point closest to zero. The bottom 
panels of Figures |3.9| and |3.10| plot the error in the numerical stationary densities 
as a function of the mesh size in logspace These numerical results confirm the 
theoretically expected rates of convergence when the mesh size is variable. 


3.7. SDEs in 2D with Additive Noise 


The examples in this part are two-dimensional SDEs of the form 

(3.7.1) dY = b{Y)dt-DU{Y)dt+^/W^dW , r(0) S D c , 

where t/ : —>■ K is a potential energy function (which will vary for each example), 

/3 is the inverse ‘temperature’, and b{x) is one of the following flow fields: 


(3.7.2) b{x) e 



= (a:i,a;2)^ 


rotational extensional shear 


nonlinear 


where 7 is a flow strength parameter. In the flow-free case (7 = 0) and in the 
examples we consider, the generator of the SDE in (3.7.1 1 is self-adjoint with respect 
to the probability density: 


v{x) = Z ^ exp(— /3C/(a:)) , Z = / exp{—pU{x))dx . 

Jn 

Moreover, the process Y is reversible m- Realizable discretizations have been de¬ 
veloped for such diffusions with additive noise in arbitrary dimensions, see [T04ll85] . 
With flow (7 7 ^ 0), however, the generator is no longer self-adjoint, the equilibrium 
probability current is nonzero, and the stationary density of this SDE is no longer 
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p(x), and in fact, is not explicitly known. In this context we will assess how well 
the second-order accurate generator Qc in (2.4.7) reproduces the spectrum of the 
infinitesimal generator of the SDE L. Since the noise is additive, the state space 
of the approximation given by Qc is a grid over C Moreover, this gener¬ 
ator can be represented as an infinite-dimensional matrix. In order to use sparse 
matrix eigensolvers to compute eigenvalues of largest real part and the leading left 
eigenfunction of Qc, we will approximate this infinite-dimensional matrix using a 
finite-dimensional matrix. For this purpose, we truncate the grid as discussed next. 


Remark 3.7.1 (Pruned Grid). Let {xij G | {i,j) G IP'} be the set of all grid 
points and consider the magnitude of the drift field evaluated at these grid points: 
{|/r(a;ij)|}. In order to truncate the infinite-dimensional matrix associated to Qc, 
grid points at which the magnitude of the drift field is above a critical value E* are 
pruned from the grid, and the jump rates to these grid points are set equal to zero. 
The price of this type of truncation is an unstructured grid that requires using an 
additional data structure to store a list of neighbors to each unpruned grid point. 
To control errors introduced in this pruning, we verify that computed quantities do 
not depend on the parameter E*. 


Remark 3.7.2 (Error Estimates). In the linear context, we compare to the ex¬ 
act analytical solution. In the nonlinear context, and in the absence of an analytical 
solution, we measure error relative to a converged numerical solution. 


3.7.1. Planar, Asymmetric Ornstein-Uhlenbeck Process. Consider (3.7.1) 
with the following two dimensional potential energy: 


(3.7.3) 


U[x) = \{x\+xl) , X = {XI,X 2 )^ G 


and the linear flows defined in (3.7.2). We will use this test problem to assess the 
accuracy of the generator Qc in representing the spectrum of L. 

In the linear context (3.7.1) with /3 = 2, the SDE (3.7.1) can be written as: 

(3.7.4) dY = CYdt + dW 


where (7 is a 2 x 2 matrix. The associated generator L is a two-dimensional, non- 
symmetric Ornstein-Uhlenbeck operator. Let Ai and A 2 denote the eigenvalues of 
the matrix C. If 3?(Ai) < 0 and 3 ?(A 2 ) < 0, then L admits a unique Gaussian 
stationary density v{x) ES See Section 4.8], and in addition, the spectrum of L in 
LP (for 1 < p < 00 ) is given by: 


(3.7.5) 


ct(L) = {niAi-I-n2A2 : ni,n2GNo}. 


The original statement and proof of (3.7.5) in arbitrary, but finite, dimension can 
be found in M. This result serves as a benchmark for the approximation. In 
the numerical tests, we select the flow rate to be 7 = 1 / 2 , which ensures that the 
eigenvalues of the matrix C have strictly negative real parts for each of the linear 
flow fields considered. 

Figure [3. ll| plots contours of the free energy of the numerical stationary density 
at 5x = 0.3 for the linear flow fields considered. At this spatial step, the error 


in the stationary densities shown is about 1%. Figure 3.12 plots the (.1 difference 
between the numerical and exact stationary densities. Figure |3.13| plots the first 
twenty eigenvalues of L as predicted by (3.7.5). Figure 3.14 plots the relative error 
between the twenty eigenvalues of largest real part of L and Qc- The eigenvalues 
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of L are obtained from (3.7.51, and the eigenvalues of Qc are obtained from using 


3.7.1 


a sparse matrix eigensolver and the truncation of Qc described in Remark 
This last figure suggests that the accuracy of Qc in representing the spectrum of L 
is, in general, first-order. 

3.7.2. Maier-Stein SDE: Locally Lipschitz Drift & Nonlinear Flow. 


The Maier-Stein system is governed by an SDE of the form (3.7.11 with fl = 


the nonlinear flow field in (|3.7.2|), /3 = 1, and 
(3.7.6) 


U{.) = 


+ 


M- 


= (xi,X2)'^ G 


2 2 

In the flow-free case (7 = 0), the Maier-Stein SDE describes the evolution of a 
self-adjoint diffusion. In this case, the governing SDE is an overdamped Langevin 
equation with a drift that is the gradient of a double well potential with two minima 
located at the points (±1,0)^. When 7 > 0, this process is no longer self-adjoint 
and the drift cannot be written as the gradient of a potential |92L 1931194) . Our 


interest in this system is to verify that (2.4.71 is convergent for SDEs of the type 
(3.7.11 with a nonlinear flow field. Eigure 3.15 | repr oduces slices of the stationary 
density, which is consistent with Figure 5 of Tm]. Contours of the free energy 
of the statio nary density at 7 = 0 and 7 = 4 are shown in Figure |3.16[ Finally, 
Figure 3.17 shows that the spectrum of L is accurately represented by Qc- 


3.7.3. Square Well SDE: Discontinuous Potential. Here, we test Qc 


(2.4.71 on a non-self-adjoint diffusion whose coefficients have internal discontinu¬ 


ities. Specifically, consider (3.7.1) with the linear flow fields in (3.7.2), /3 = 1, and 


a regularized square well potential: 


(3.7.7) 


UM = ,a„h - _ 7 , 


— tanh 


^max(a;i, 3:2) — ^2 


= (X1,X2)'^ G fl = [-1,1] X [-1,1] 


where di and d 2 specify the location of the discontinuity and e is a smoothness 
parameter. We assume that the solution is periodic at the boundaries of Q and we 
discretize the domain so that all of the grid points are in the interior of H. Since 


the potential energy is nearly discontinuous, we replace the jump rates in (2.5.2) 
with the following rates which avoid differentiating the potential energy 


(3.7.8) 


q{x,y) = ^exp (-^{U{y) -U{x) - b{x)'^ {y - x)) 


To be clear, x, y in (|3.7.8|) are 2-vectors. Contours of the numerical stationary 
density are plotted in Figure [3.18| when the stationary density is known. The spatial 
grid size used is h = 0 . 1 , which is selected so that the numerical densities represent 
the true stationary densities to within 1% accuracy in the £i norm. Figure 3.19 


shows a similar plot in the flow cases for which the stationary densities are unknown. 
In this case we checked that at the spatial grid size of h = 0.1 the stationary 
density is well-resolved. To accurately represent the probability flux through the 
square well, we construct the coarsest cell-centered mesh to have borderlines that 
coincide with the discontinuities of the square well. Subsequent refinements involve 
halving the size of each cell, so that the discontinuities always lie on a cell edge. 
Figure |3.20| illustrates the £2 convergence of the numerical estimate to the first 
twenty eigenvalues of the continuous operator L. As a side note, one may be able 
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to construct a more accurate generator in this context by extending the spatial 
discretization introduced in HI] to non-self-adjoint diffusions. 


3.7.4. Worm-Like Chain SDE: Singular Drift. Consider (3.7.1) with the 


linear flow fields in (3.7.2), /3 = 1, and the potential energy given by: 


(3.7.9) U{x) = 


1 


1 - Ixil 


- Ixil -I- 2x1 


X = ixi,x2)'^ en = (- 1 , 1 ) 


This system is an example of an SDE with a boundary and with a drift vector field 
that has singularities at the boundaries of D: (± 1 , 3 : 2 )^ for all X 2 € K. Because 
of this singularity, standard integrators for ( 3.7.1| ) may diverge in finite time even 
in the flow-free case. Contour lines of the numerical stationary density i'(x) are 


plotted in Figure 3.21 and 3.22 when the stationary density is known and unknown, 
respectively, and for a flow rate of 7 = 5. We choose a high flow rate to make the 
numerical test more rigor ous. The accuracy of the spectrum of the generator Qc 


is assessed in Figure 


3.23 


This figure suggests that is second-order accurate in 


representing the spectrum of the generator L. 


3.8. Adaptive Mesh Refinement in 2D 

For the subsequent planar SDE problems with multiplicative noise, we use a 
variable step size grid with the property that the generator in (2.4.7) evaluated 
on this grid is realizable. Recall, if the diffusion matrix is diagonally dominant, 
a grid with equal spacing in the horizontal and vertical directions suffices. If the 
diffusion matrix is constant, but not diagonally dominant, the ratio of the spacing in 
the horizontal and vertical directions can be chosen so that the grid has the desired 
property, since a 2 x 2 , symmetric, and positive-definite matrix is always weakly 
diagonally dominant [8|. In general, if the diffusion matrix is neither diagonally 
dominant nor constant, it is not straightforward to construct a grid such that Qc 
is realizable on this grid. However, in the specific context of the log-normal and 
Lotka-Volterra examples considered in this paper, the following lemma states that 
a grid with the desired property exists. The key idea in this construction is to make 
the grid evenly spaced in logspace. 

To state this lemma, define: 


diag(a:) = 


xi 0 
0 X2 


for any x = {xi,X 2 )"^ G 


Lemma 3.8.1. Consider the generator Qc in (2.4.7) with domain Ki, an arbi¬ 


trary drift field, and a diffusion matrix of the form: 


(3.8.1) 

Assume that 

(3.8.2) ^ 0 


M{x) = diag(x) 


Mil J^12 

M12 j^22 


diag(a:) . 


det 


Mil J^12 

M 12 ^22 


> 0 , Mil J^22 ^ Q ^ 


In this context, Qc in (2.4.7) is realizable on a gridded state space. 
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Proof. To prove this lemma we construct a grid such that Qc in (2.4.7) is 


realizable on this grid. Let S = {{xi,yj)} C be the desired grid. We define S 
to be the image of an evenly spaced grid in logspace: 


S = e 


= Ci+1 = Vj+1 - Vj} 


where and Srj set the grid size in the horizontal and vertical directions in logspace, 
respectively. We choose the spacing in logspace so that S has the desired property 
in physical space. Specifically, map every point G S' to a point in (xt, yj) G S 

using the transformation: Xi = exp(^i) and yj = exp{r]j). Suppose that 


= ae and Srj = e 


where e and a are positive parameters. We choose these parameters so that Qc 
evaluated on S is realizable. Note that the spacing between neighboring horizontal 
and vertical grid points in physical space satisfy: 


(3.8.3) 


Sxf = Xi+i - Xi = (exp(ae) - l)xi 
Sx~ = Xi — Xi-i = (1 — exp(—ae))a:i 
5x, = t = sinh(ae)a;, 


and 


(3.8.4) 


<5?// = y]+i - Vj = (exp(e) - l)yj 
^V7 = Vj - Vj-i = (1 - exp(-e))yj 


= — 


Sy^ 

—^ = sinh(e)i/j 


Referring to (2.4.7), a necessary and sufficient condition for realizability of Qc 


is that all of the following inequalities hold: 


M 


11 


SxiSxf 


Ml] V 0 
6xfSyl 


Mil AO 

hJ > Q 


5x+Sy- 


Ml] Ml] V 0 Ml] A 0 

Sxi5x~ 6x~6y~ 6x~6yj' 

Ml] Ml] V 0 Ml] A 0 

SyjSy]' SxfSyl' Sx~dy]' 

Ml] Ml] VO Ml] AO 

SyjSyJ Sx~dy~ Sxf 6y~ 


(3.8.5) 


















40 


3. EXAMPLES & APPLICATIONS 


Let US express these inequalities in terms of a. and e. Substitute ( |3.8.1[ ) and the 

d simplif; 

AO 


grid spacings (3.8.3) and (3.8.4) into (3.8.5) and simplify to obtain: 


(3.8.6) 


VO 

sinh(ae) exp(e) — 1 1 — exp(—e) 

V 0 A 0 

sinh(ae) 1 — exp(—e) exp(e) — 1 

^22 ^12 Q ^12 Q 


sinh(e) exp(ae) — 1 1 — exp(—ae) 

M22 V 0 A 0 

sinh(e) 1 — exp(—ae) exp(ae) — 1 

Applying the elementary inequality 

exp(a) — 1 > 1 — exp(—a) V a S 1 


> 0 

> 0 

> 0 

> 0 


to (3.8.6) yields the following sufficient condition for realizability 


(3.8.7) 


|Mi 


> 


sinh(ae) 

1 — exp(—e) 


and 


1 — exp(—ae) \M 


sinh(e) 


> 


M 22 


The hypothesis in (3.8.2) implies that 


> 0 , M22 > 0 and 


M 


11 


\M 


121 


\M 


121 


M 22 


> 0 . 


Hence, we can choose a to satisfy: 


\M 


121 


|M12| 

> “ > ^ ■ 


This choice of a implies (3.8.7) holds if e is small enough. Indeed, a straightforward 
application of I’Hopital’s Monotonicity Rule implies that 

1 — exp(—ae) 


(3.8.8) 


—, . sinh(ae) 

= 1 - exp(-e) 


> a > 


= <('(e) 


sinh(e) 

for any e > 0. Specifically, to prove the upper bound, note that: 

lim 0(s) = a . 

s —^0 


Moreover, I’Hopital’s Monotonicity Rule implies that the function (/)(s) is increasing 
since the function: 

s I—>■ a cosh(as) exp(s) 

is increasing with s. Thus, 0(s) is increasing. Similarly, ^(s) is decreasing since the 
function: 

s I—>■ aexp(—as) sech(s) 

is decreasing with s, and lims_,.o (/'(s) = a. □ 
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3.9. Log-normal Process in 2D with Multiplicative Noise 

We derive a log-normal process by transforming a two-dimensional Ornstein- 
Uhlenbeck process. Thus, we obtain a log-normal process whose transition and 
equilibrium probability distributions are explicitly known and can be used as bench¬ 


marks to validate the approximation given by Qc in (2.4.7). First, let us recall some 


basic facts about Ornstein-Uhlenbeck processes. Consider a two-dimensional, non- 
symmetric Ornstein-Uhlenbeck process 

(3.9.1) dY = AYdt + V2BdW , Y{0)=xeR^ . 

Let M be the diffusion matrix defined as: 


M = BB^ = 




M22 


Recall that the transition probability of Y(t) is a two-dimensional Gaussian distri¬ 
bution with mean vector and covariance matrix S(t) given by: 

^{xA) = ex.p{tA)x , 


E(t) = 2 / exp(sA)M exp(sA^)ds 

Jo 


See, e.g., [en Section 4.8]. To numerically compute E(t) for any t > 0, observe 
that the covariance matrix S(t) satisfies the Lyapunov equation: 

T,{t)A'^ + A^{t) = 2 y exp{sA)M exp(sA^)ds^ , 

= 2 {^exp{tA)M{exp{tA))^ — M) . 

The MATLAB command lyap can be used to solve equations of this type. Specif¬ 
ically, given the matrices A and B as inputs, lyap solves for the unknown matrix 
X which satisfies the equation: 

AX + XA^ +B = Q . 

Defining the inputs to lyap as: 

A = A , B = —2 (exp(tA)M(exp(tA))^ — M) , 

the function will output the desired covariance matrix. (A sufficient condition for 
a unique solution to this Lyapunov equation is that all of the eigenvalues of the 


matrix A have negative real parts.) If the eigenvalues of the matrix A in (3.9.1) 
all have strictly negative real part, passing to the limit in the above expressions as 
t ^ oo implies that the equilibrium distribution is Gaussian with zero mean vector 
and a covariance matrix Sqo that satisfies the Lyapunov equation: 

(3.9.2) SooA^-H AEoo =-2M . 

Gonsider the transformation of Y{t) = {Yi{t),Y 2 {t))'^ given by 

A = (exp(yi),exp(y2))^ ■ 

Under this transformation, the linear SDE ( 3.9.1| ) becomes: 

(3.9.3) dX = fi{X)dt + V2a{X)dW , A(0) S 
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where the drift and diffusion fields are: 

4 log(a:) 


fx(x) = 

M(x) = (era 

= diag(a;i,a:2) 


■^)(x) 


diag(a;, y)A 




12 


M 

M 22 


log(a; 2 ) 


diag(xi,a:2) = 


M^^xf 

M^^XiX2 


M^^XiX2 

M^^x^ 


By change of variables, the (log-normal) probability and stationary densities of the 
transformed process X(t) are given by: 

Pt(x) = pt(log(xi),log(x 2 ))(xiX 2 )^^) 

vix) = p(log(a;i),log(a;2))(aiia;2)"^ J 


(3.9.4) 


= (xi,X2y 


where ptix) and ie(x) are the (normal) probability and stationary densities of the 
Ornstein-Uhlenbeck process Y(t), respectively. 


Figure 


3.24 


illustrates stationary density accuracy of Qc in (2.4.7) using the 


variable step size grid derived in the preceding Lemma |3.8.1| 

3.10. Lotka-Volterra Process in 2D with Multiplicative Noise 

Consider (2.1.2) with 12 = Ki and set: 


p(x) = 


kixi — xiX2 — Jixl 
-k2X2 + XiX2 — 723^2 


and (aa'^)(x) = diag(a:) 


Mil J^12 

M12 j^22 


diag(x) 


This SDE is a model of interaction between a population of Xi prey and X 2 predators 
that includes environmental fluctuations. Recall that when 71 = 72 = 0, the noise- 
free system reduces to an ODE with a fixed point at (fc 2 ,Ai) and a first integral 
given by the following function: V(x) = log(xi) -I- ki log(x 2 ) — (xi -I- X 2 ). The 
diffusion field of this SDE satisfies the hypotheses in Lem ma|3.8.1| Thus, we will 
be able to construct a grid with the property that Qc in (2.4.7) evaluated on this 
grid is realizable. 

We use some known theoretical results on this SDE problem to benchmark 
the generator Qc mg. For this purpose, it helps to define the parameters: ci = 
ki — 5 Mil C 2 = ^2 -f jM^^. According to Theorem 1 of [122) . for every 
t > 0 the probability distribution of the SDE is absolutely continuous with respect 
to Lebesgue measure on D. This property implies that the solution is strictly 
positive for all t > 0. This positivity, however, does not prevent the probability 
distribution of the solution from converging in the long time limit to a measure 
supported only on the boundaries of D, in which case one or both of the populations 
may become extinct. In fact, the finite-time probability distribution of the SDE 
solution converges to an invariant probability measure (IM) with the following 
support properties: 


If Cl > 71 C 2 : then the IM is supported on 
If 71 C 2 > Cl > 0: then the IM is supported on 
xi given by: 


(3.10.1) 


Mxi) = z 


- 1 , 


exp 


X {0} with marginal density of 

2 xi \ 

Mil ) 
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where Z is the following normalization constant: 

/ 2ci 

^ J [m^ 

If Cl < 0: then the IM is atomic at the origin. 




3.25 


Note that if 2ci/M^^ < 1 (or fci < M^^), then the marg inal stationary density 

use Qc in (2.4.7) with 


and 


3.26 


in (3.10.1) is unbounded at x = 0. Figures 
the variable step size grid described in Lemma 3.8.1| to numerically verify these 
theoretical results. 


3.11. Colloidal Cluster in 39D with Multiplicative Noise 

Consider a Brownian Dynamics (BD) simulation of a cluster consisting of 13 
particles governed by the following SDE: 

(3.11.1) 

dV = -M(V)D£(Y)dt + dW M(Y)dt + y/2^a{Y)dW , r(0) G K" , 
tT(x)cr(a;)^ = M{x) , 


where £(x) is the potential energy of the system, M{x) is the mobility matrix (or 
what we have been calling the diffusion matrix), and IF is a n-dimensional Brownian 
motion. To be clear, the total dimension of this system is n = 3 x 13 = 39. This BD 
model of a colloidal gel demonstrates that hydrodynamic interactions can influence 
the collapse dynamics of the cluster |35) as verified below. 

The potential energy £{x) involves short-range, pair-wise interactions between 
particles. It features a very steep hard core with a short-range attractive tail. 
When hydrodynamic interactions are accounted for, the mobility matrix depends 
on the state of the system and its entries involve long-range interactions that decay 
like 1/r with distance r between particle pairs. The precise form of £ and M 
are derived below. In the situations we consider, the divergence of the mobility 
matrix appearing in (3.11.1) vanishes. Additionally, this diffusion is self-adjoint 


with respect to the density: 


(3.11.2) 


iy{x) = exp{—/3£{x)) . 


For background on the numerical treatment of self-adjoint diffusions with multi¬ 
plicative noise see m 

For this BD problem, numerical stability and accuracy impose demanding time 
step size requirements for time integrators due to the presence of singular coefficients 
in the drift of the governing SDE. Another issue with this SDE problem is that the 
hydrodynamic interactions between particles induce long-range interactions, which 
leads to an SDE with multiplicative noise. 

Potential Energy. Consider a collection of N identical spherical particles with 
radius a. Let q = (qi,... ,qN) G denote the position of all N particles where 
qi G for i = 1,..., N. The distance between the ith and jth particle is the usual 
Euclidean one: 


(3.11.3) 


d{qi,qj) = h -qj \ ■ 
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In terms of this distance, the total energy is a sum of pair-wise potentials over all 
particle pairs: 


(3.11.4) 


£{q) = 


E 

l<i<n 

j>i 


U{d[qi,qj)) , 


Following the description in 
U(r) is given by: 

(3.11.5) U{r) = Usc{r) + UAo{r) . 

which is a sum of a soft-core potential: 

24 


, we assume that the pair-wise potential energy 


(3.11.6) 


C^sc(’') = esc ( — 
r 


and an attractive potential that is defined in terms of the inter-particle force Fao (^) 
it exerts as follows: 

(3.11.7) 

{cAo{Dl-Dl) iir<Di 

UAo{f) = J FAo{s)ds , PAoif) = < cao{D2 — r^) ii Di < r < D 2 

I 0 otherwise 


Here, we have introduced the parameters esci CAO, Di, and £> 2 . Figure [3.27| plots 
graphs of C/(r) and F{r) for the parameter values given in Table[^ Note that U{r) 
induces a short-range force. Further, because of the 1/r^^-singularity in Uscif) 
in ( 3.11.6| ) at r = 0, this force makes the resulting SDE stiff. Since the energy is 
short-range, it is not uniformly coercive. Thus, the invariant density in (3.11.2) 
is not integrable on K" and the system is not ergodic. In the numerical tests, we 
observe that - despite this non-ergodicity - for most sample paths, the cluster does 
not evaporate over the time-span of simulation. 

Mobility Matrix. Consider again a collection of N identical spherical particles 
with radius a. Here we describe a commonly used, implicit model for the effect of 
the solvent known as the Rotne-Pragner-Yamakawa (RPY) approximation. This 
model results in solvent-mediated interactions between the N spherical particles 
through the mobility matrix M{q). Let /ax 3 be the 3x3 identity matrix, i?hydro be 
the hydrodynamic radius of each bead, and rjs be the solvent viscosity. The RPY 
mobility matrix is given by: 


^/3x3 , if*=j 

^Rpviqi — Qj) , otherwise 




^1,N 


(3.11.8) M{q) = 



J 


^N,l ■ 

■ £iN,N 

1 


Here, we have introduced the 3x3 matrix flupy^x) defined as: 
£Irpy{x) = ^ ^Ci(a;)/3x3 + ^2(a:)|^ O 


for all q G 
(3.11.9) 

where C'i(x) and (72 (x) are the following scalar-valued functions: 

I ( "ixt"") ”*" 5 1^1 ^ 2i?hydro 


Ci(x) = 


1 - 
^ 32 


\ ^hydro / ’ 


otherwise 
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and, 


C2ix) = 


( -Rhydr< 

“kT 


y Rhydro J ’ 


/ Rhydro \ 

V l-l ) 


if |X| > 

otherwise 


The quantity 1/C is the mobility constant produced by a single bead translating in 
an unbounded solvent at a constant velocity: C = 67r77si?hydro- The approximation 
(3.11.8) preserves the physical property that the mobility matrix is positive definite, 
satisfies {div M){q) = 0 for all q S , and is exact up to 0((i?hydro/Tij )^) where 
rij is the distance between distinct particles i and j. To read more about the RPY 
approximation see [m]. 

Description of Numerical Test. We are now in position to repeat the simulation 
using the generator Qc in (2.5.2). Consider 13 particles initially placed 


on the vertices and center of an icosahedron of edge length 8.08, which corresponds 
to a radius of gyration of 7.68. Note that this edge length is slightly greater than 
the equilibrium bond length, which is 6.4. Thus, without noise, the cluster would 
simply relax to an icosahedron with an edge length equal to the equilibrium bond 
length. To prevent exploding trajectories that can occur if any two particles get 
too close, we set the spatial grid size to be 10% of the equilibrium radius or 0.32. 
This choice ensures that the numerical approximation to the repulsive force cannot 
cause any particle to take a single spatial step that is larger than the cutoff radius 
of the short-range interaction. Thus, the approximation to the repulsive force 
between every particle pair is sufficiently resolved in space to prevent explosion. 
Figure and [3.29| illustrate the results produced by an SSA integrator induced 
by the generator Qc in ( 2.5.2[ ) operated at this spatial step size. Other simulation 
parameters are given in Table An important physical parameter is the Brownian 
time-scale defined as: ts = prjsa^. In order to compare to the data from | 25 ) . the 
time span of simulation is long and set equal to T = dOOts = 1065. We emphasize 
that the spatial grid size with and without hydrodynamic interactions is set equal 
to 10%a or 5x = 0.32. Finally, Figure 3.30 validates Proposition 4.3.11 which 


quantifies the average number of computational steps of the SSA integrator as a 
function of the spatial step size. 
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Parameter 

Description 

Value (s) 

Physical Parameters 

N 

# of particles 

13 

Vs 

solvent viscosity 

1 

r' 

temperature factor 

12.3 

a 

particle radius 

3.2 

^hydro 

hydrodynamic radius 

{«, 0} 

esc 

parameter in Use 

10 

Di 

parameter in Uao 

2.245a 

D2 

parameter in Uao 

2.694a 

CAO 

parameter in Uao 

58.5/0^ 

Numerical Parameters 

5x 

spatial step size 

10% a 

T 

time-span of simulation 

1065 


Table 1. Colloidal Collapse Simulation Parameters. Physical 
parameter values are taken from: [3511^ . 
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spatial step size 




spatial step size 


Figure 3.2. Cubic Oscillator. From top: accuracy of the gen¬ 
erators Qu in ( |2.3.6[ ) and Qc in ( |2.3.9| with respect to the stationary 
density, mean first passage time, and committor function of the SDE. A 
formula for the true solutions is explicitly available in this example and 
serves as a benchmark for comparison. For the numerical solution, we 
use the formulas given in Sections |6.1| |6.4| and |6.3| In all cases we see 
that Q„ (resp. Qc) is first (resp. second) order accurate in representing 
these quantities. The reference solutions are graphed in the insets. 
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3. EXAMPLES & APPLICATIONS 


Sample Paths 



Figure 3.3. Cubic Oscillator Simulation. This figure plots paths 
produced by the SSA induced by (light grey), Qc (grey) and Qe 
(black). In all cases the meshes are eqnidistant with 5x = 0.25 and the 
initial condition is marked by the black dot at t = 0. The paths are 
driven by the same sequence of nniform random variables. We stress 
that the SSA produced by Qe uses the exact mean holding time of the 
SDE solution, which far enough from the origin is well approximated 
by ( |3.3.4[ ). The figure in the inset shows how the SSA time evolves 
with the SSA step index. Far away from the origin, the dynamics is 
dominated by the deterministic drift, which causes the walkers to move 
toward the origin and then become localized around the origin. The SSA 
associated to Qc moves to the origin too fast because its mean holding 
time tends to zero too fast, as (3.3.61 predicts. In contrast, the mean 
holding times for Qu and Qe are in close agreement far from the origin 
(where the asymptotic analysis in (3.3 is valid), but closer to the origin, 
this agreement is not apparent, as expected. 
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Figure 3.4. Log-normal Process in ID. This figure confirms that 
the realizable discretizations given in (3.1.11 and (3.1.21 are accurate 
with respect to the mean first passage time (bottom panel) and exit 
probability or committor function (top panel) computed on the interval 
[1/2, 5]. The reference solution is shown in the inset. An evenly spaced 
mesh is used to produce this numerical result, with spatial step sizes as 
indicated in the figures. The figure confirms that (resp. Qc) is first 
(resp. second) order accurate. 
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spatial step size 


Figure 3.5. Log-normal Process in ID. The top panel plots the 
free energy of the stationary density of the SDE. The bottom panel 
plots the accuracy of Qu and Qc in representing this stationary density 
as a function of the spatial step size in log-space. The benchmark 
solution is the stationary density (3.5.11 evaluated on the variable step 


size grid described in i 3.4 Note that the central scheme Qc appears 


to be second-order accurate in representing the stationary probability 
distribution of the SDE, which numerically supports Theorem |4. 5. 2[ 
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Figure 3.6. Log-Normal Process in ID. The top panel plots 
realizations of the generators Qu, Qc, and Qe produced using the SSA 
integrator with the initial condition at X(0) = 20 marked by the black 
dot at t = 0. Recall that Qe reproduces the true mean holding time 
of the process, see (2.3.151 for its definition. The state space of the 


processes is the variable step size grid described in §3.4| with a grid size 
in log-space of 5^ = 0.25. The time interval of simulation is [0,100]. 
The mean SSA time for all three methods is approximately 0.03. The 
bottom panel plots the SSA time as a function of the computational 
step index. This figure illustrates that the SSA integrators never exit 
the domain of definition of the coefficients of the SDE. 
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3. EXAMPLES & APPLICATIONS 




Figure 3.7. Log-Normal Process in ID. The top (resp. bottom) 
panel plots the accuracy of the generators Qu and Qc with respect to 
Ea:(F(t)^) as a function of spatial step size (resp. mean holding time) 
with initial condition y(0) = 2 and terminal time t = 1. Note that the 
central scheme Qc appe ars to be second-order accurate, which numeri¬ 
cally supports Theorem 4.5.1 This figure was produced using 10® SSA 
realizations. 
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Figure 3.8. Log-Normal Process in ID. These statistics are 
produced using the SSA induced by the generators Qu in (3.1.11 and Qc 
in ( |3.1.2[ ) as indicated in the figure Legend. The top panel of this figure 
plots the mean number of computational steps of the SSA integrator as 
a function of the spatial step size. The bottom panel of this figure plots 
the mean holding time of the SSA integrator as a function of the spatial 
step size. The statistics are obtained by taking a sample average over 
100 realizations with initial condition at A'(O) = 1 and a time interval 


of simulation of [0,10]. This figure validates Proposition 4.3.11 
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Figure 3.9. Cox-Ingersoll-Ross Simulation. The top panel plots 
the free energy of the stationary density of the SDE, when the boundary 
is natural. The bottom panel plots the accuracy of and Qc in repre¬ 
senting this stationary density as a function of the spatial step size in 
log-space. The benchmark solution is the cell averaged stationary den¬ 
sity (3.6.31 using the variable step size grid described in [ 3.4 Note that 


the central scheme Qc is second-order accurate in representing the sta¬ 
tionary probability distribution of the SDE, which numerically supports 
Theorem 14.5.21 
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Figure 3.10. Cox-Ingersoll-Ross Simulation. The top panel 
plots the free energy of the stationary density of the SDE, when the 
boundary is regular. The bottom panel plots the accuracy of Qu and 
Qc in representing this stationary density as a function of the spatial 
step size parameter 5^ in log-space. The benchmark solution is the 
cell averaged stationary density ( |3.6.3[ ) using the variable step size grid 
described in i3.4 Note that the central scheme Qc is second-order 


accurate in representing the stationary probability distribution of the 
SDE, which numerically supports Theorem |4. 5. 2| 
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Figure 3.11. Planar, Asymmetric Ornstein-Uhlenbeck Pro¬ 
cess. These figures plot contours of the free energy of the numerical 
stationary density associated to Qc in the flow-free (a), rotational (b), 
extensional (c), and shear flow (d) cases. The grid size is selected to be 
h = 0.3, which yields an approximation to the true stationary density 
that is within 1% error in the £i norm. The circular contour lines of 
the stationary densities in the (a) flow-free and (b) rotational cases are 
identical. In the (c) extensional and (d) shear flow cases, the contours 
are elliptical. The flow strength parameter is set at: 7 = 1/2. 
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spatial grid size 


Figure 3.12. Planar, Asymmetric Ornstein-Uhlenbeck Pro¬ 
cess. These figures plot the error in the estimate of the true sta¬ 
tionary density as a function of grid size for in the flow-free (FO), 
rotational (FI), extensional (F2), and shear (F3) flow cases, as indicated 
in the figure legend. The figure shows that this approximation is second- 
order accurate in estimating the stationary density. The flow strength 
parameter is set at: 7 = 1/2. This numerical result is consistent with 
Theorem 14.5.21 
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Real Eigenvalues 


Complex Eigenvalues 
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10 15 

n 

(d) shear 


20 


Figure 3.13. Planar, Asymmetric Ornstein-Uhlenbeck Pro¬ 
cess. These figures plot the twenty eigenvalues of largest nontrivial real 


part of the operator L as predicted by the formula given in (3.7.51. In 


(a), (c) and (d) the formula implies that the eigenvalues are strictly 
real, while in the rotational flow case shown in (b) the eigenvalues are 
complex. The flow strength parameter is set at: 7 = 1/2. These eigen¬ 
values serve as a benchmark solution to validate the accuracy of Qc in 
representing the spectrum of L. 
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Accuracy in a[L) 



Figure 3.14. Planar, Asymmetric Ornstein-Uhlenbeck Pro¬ 
cess. These figures plot the relative error of the numerical estimate 
produced by Qc of the twenty eigenvalues of largest nontrivial real part 
of the operator L in the flow-free (FO), rotational (FI), extensional (F2), 
and shear (F3) flow cases, as indicated in the figure legend. The figures 
suggest that the accuracy of Qc in representing these eigenvalues is 0{h) 
in general, even though it is 0{h^) in the flow-free (FO), rotational (FI), 
and extensional (F2) flow cases. The flow strength parameter is set at: 
7 = 1/2. 
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Figure 3.15. Maier-Stein: Slices of Stationary Density. These 
figures plot the intersection of the stationary density of the Maier-Stein 
system with the planes indicated in the legends with parameters: = 2 
and 7 = 4.67. These graphs agree with those given in Figure 5 of [ml. 
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(a) flow-free (b) nonlinear flow 


Figure 3.16. Maier-Stein: Stationary Densities. These figures 
plot contours of the free energy of the stationary density for the gener¬ 
ator Qc with grid size h — 0.05. The test problem is a particle in the 
double-well potential given in (3.7.61 and the nonlinear flow field given 
in (3.7.21. The flow strength parameters are set at: 7 = 0 (left panel) 
and 7 = 4 (right panel). 



spatial grid size 


Figure 3.17. Maier-Stein: Accuracy in o-{L). This figure illus¬ 
trates the convergence in the £2 norm of the numerical estimate produced 
by Qc of the twenty eigenvalues of largest nontrivial real part of the op¬ 
erator L for the values of the flow strength parameter 7 indicated in the 
legend. 
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8 


6 


4 


(a) flow-free 


(b) extensional 


Figure 3.18. Square Well: Stationary Densities (known). 

These figures plot contours of the free energy of the numerical station¬ 
ary density for a particle in the potential given by (3.7.71 with j3 = 1, 
di = 0.4, d 2 = —0.4, e = 0.001, and the flow-free (left panel) and ex¬ 
tensional flow (right panel) cases with 7 = 2. For these flows a formula 
for the stationary density i'{x) is available, and the numerical densities 
shown are within 1% accurate in the ^inorm. 
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(a) rotational 


(b) shear 


Figure 3.19. Square Well: Stationary Densities (unknown). 

These figures plot contours of the free energy of the numerical station¬ 
ary density for a particle in the potential given by (3.7.71 with j3 = 1, 
di = 0.4, d 2 = —0.4, e = 0.001, and the shear (left panel) and rota¬ 
tional flow (right panel) cases with 7 = 2. For these flows a formula for 
the stationary density v{x) is not available, but we verified that these 
numerical solutions are well-resolved. 
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Figure 3.20. Square Well: Accuracy in o{L). This figure plots 
the relative error in the £2 norm of the numerical estimate produced 
by Qc of the twenty eigenvalues of largest nontrivial real part of the 
operator L. The parameters here are set at: 7 = 1 and j3 = 1. The 
figure suggests that the relative error of Qc in representing the spectrum 
of L is about second-order in the spatial step size. The slight divergence 
at small spatial step size is attributed to errors in the sparse matrix 
eigensolver used. 








Figure 3.21. Worm-Like Chain: Stationary Densities 

(Known). These figures plot contours of the stationary density for a 
particle in the potential given by (3.7.91, and from left, flow-free and the 
extensional flow given in (3.7.21 with 7 = 5 and 13 = 1. For these flows 
a formula for the stationary density i'{x) is available, and the numerical 
densities shown are within 1% accurate in the ^inorm. 
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Figure 3.22. Worm-Like Chain: Stationary Densities (Un¬ 
known) . These figures plot contours of the free energy of the numerical 
stationary density for a particle in the potential given by (3.7.91, and 
from left, rotational and shear flows given in (3.7.21 with 7 = 5 and 
P = l. 
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spatial grid size 


Figure 3.23. Worm-Like Chain: Accuracy in t(L). This figure 
plots the relative error in the I 2 norm of the numerical estimate produced 
by Qc of the twenty eigenvalues of largest nontrivial real part of the 
operator L. Parameters are set at 7 = 1 and /3 = 1. This figure suggests 
that Qc is about second-order accurate in representing the spectrum of 
the generator L. 
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Figure 3.24. Log-normal Process in 2D. This figure illustrates 
stationary density accuracy of in (2.4.71. The grid nsed is evenly 
spaced in logspace as shown in the upper left inset. Contonrs of the 
stationary density are plotted in the lower left inset. The spatial step 
size shown is with respect to the spacing of the grid in logspace. In this 
example the benchmark solution is the true stationary density given by 
the formula in (3.9.41. See Lemma 3.8.1 for more details about the grid. 
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Figure 3.25. Lotka-Volterra (LV) Process. This figure plots 
contours of the numerical stationary density for the LV process pro¬ 
duced by the generator Qc in (2.4.71 with the variable step size grid 
described in Lemma [3.8.1| The parameters are selected so that the true 
stationary distribution is supported on R+. This figure is in agreement 
with Theorem 1 of [T22]. 
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Stationary Density Plot, Ci =0.74, C 2 =1.02 



Figure 3.26. Lotka-Volterra (LV) Process. Here we numerically 
verify Theorem 1 of [T22] using the generator Qc in ( |2.4.7| with the 
variable step size grid described in Lemma The figure plots the 

density of the stationary distribution for the LV process in a parameter 
regime where it is supported on R+ x {0}. In this case in the long 
time limit, the ‘predator population’ (or X 2 degree of freedom) becomes 
extinct, and an explicit formula for the stationary density in the ‘prey 
population’ (or xi degree of freedom) is given in (3.10.11, which we use 
to benchmark the numerically produced solution. 
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Interparticle Force & Potential 



Figure 3.27. Graph of pair-wise potential energy and force. 

Graph of the inter-particle energy U{r)/cAO given in (3.11.51 and corre¬ 
sponding inter-particle force F{x)/cao = —U'if )/cao as a function of 
the inter-particle distance r with cao = 1.7853 and a = 3.2. The range 
of interaction is specified by the parameter D 2 , which here is set equal 
to D 2 = 2.694a « 8.6. The energy is minimized at twice the particle 
radius a with minimum value: U{2a) ~ —51.82. Up to a factor cao, this 
graph agrees with Fig. SI in the Supplementary Material accompanying 










72 


3. EXAMPLES & APPLICATIONS 



step index 


Figure 3.28. Brownian Dynamics Simulation in 39D. This 
figure shows sample paths generated by an SSA integrator of the col¬ 
lapse of a colloidal cluster. The top panel shows a sample trajectory 
with and without hydrodynamic interactions (HI). The bottom panel 
shows the SSA time vs step index for these sample trajectories. These 
sample paths are produced using the SSA induced by the generator Qc 
in (2.5.21 operated at a spatial step size set equal to 10%a, where 2a is 
the equilibrium bond length from Figure [3.27| 
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Figure 3.29. Brownian Dynamics Simulation in 39D. This 
figure shows the mean radius of gyration of the colloidal cluster as a 
function of time over the interval [0, T] computed using an SSA inte¬ 
grator (solid) and a benchmark (dashed) from [35j with and without 
hydrodynamic interactions (HI). This numerical result is produced us¬ 
ing the SSA induced by the generator Qc in (2.5.21 operated at a spatial 
step size set equal to 10%a, where 2a is the equilibrium bond length 
from Figure [T^TI 
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Figure 3.30. Complexity of SSA in Brownian Dynamics Sim¬ 
ulation in 39D. The top panel of this figure plots the mean number of 
computational steps of the SSA integrator as a function of the spatial 
step size. The bottom panel of this figure plots the mean holding time 
of the SSA integrator as a function of the spatial step size. The figure 
validates Proposition |4.3.11| This numerical result is produced using 
the SSA induced by the generator Qc in (2.5.21. 


















CHAPTER 4 


Analysis on Gridded State Spaces 


In this chapter we analyze the stability, complexity, and accuracy of the real¬ 
izable discretization Qc in (2.5.2) on a gridded state space. (Refer to (2.2 for the 
definition of a gridded state space.) For the stability analysis, we assume that the 
drift field is locally Lipschitz continuous and weakly dissipative. To show that the 
approximation is second-order accurate, we further assume that the drift field is of 
class C"*. In order to ensure that the state space is gridded, we assume that the 
noise is additive and isotropic. These assumptions are sufficient to prove that the 
approximation has a stochastic Lyapunov function, which is then used to quantify 
the complexity and accuracy of the approximation. 


4.1. Assumptions 


Here we describe the main objects of the theory in this chapter: 




generators of the processes, given in (4.1.4) and (4.1.8); 

stochastic equations governing the dynamics of sample paths, given in 

diTsj ) and ( |4.1.12[ ); and, 

Kolmogorov equations governing the evolution of conditional expectations 
of an observable, given in (4.1.6) and (4.1.14). 


4.1.1. Some Notation. We shall use the following notation. Let | • | denote 
the 2-norm of an n-vector, let | • |p denote the p-norm of an n-vector where p >1 and 
let 11 • II denote the standard operator norm. For a real-valued function / G C"’(K") 
where r G N, let £>’'/ denote the rth derivative of /. Let denote the Banach 

space of all Borel measurable and bounded functions from K" to K endowed with 
the supremum norm: 

(4.1.1) ll/lloo = sup \f{x)\ , / G Hb(K") . 

xGR" 

We denote by Cb{R^) the subspace of Hb(M") consisting of all uniformly continuous 
and bounded functions. This subspace with the supremum norm: 

(4.1.2) ll/lloo = sup |/(a:)| , / G ^(R") , 

a;GK" 

is a Banach space. For A: G N, we denote by (^^(R”) the subspace of all fc-times 
differentiable functions whose derivatives up to fcth order are uniformly continuous 
and bounded in the supremum norm. This subspace is also a Banach space when 
endowed with the norm: 

k 

(4.1.3) ll/IU = ll/lloo+ ^ sup ||Z17(*)II , /GC,^(R7. 

. , xGK" 
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4.1.2. SDE. Consider an infinitesimal generator L, which maps a twice dif¬ 
ferentiable test function / : K” —)■ K to another function Lf : K” —)■ K defined 
as: 


Lf{x) = trace(I?/(a;)/r(a:)^ -I- D‘^f{x)) 


(4.1.4) 

where /r : K" —>■ K" is the drift field. The SDE associated to L has an unbounded do¬ 
main K" and isotropic (or direction-independent), additive (or state-independent) 
noise: 


(4.1.5) 


dY = fi(Y)dt + V2dW , r(0) = a; e 


Under certain assumptions on the drift field (e.g. those stated in Assumption |4.1[ ), 
this SDE uniquely defines in a pathwise sense a Markov process Y(t) for any t > 0. 
Basic properties of solutions to (4.1.5|) are reviewed in (4.3 


For any < > 0 and x G K", let Pt and \it,x resp. denote the Markov semigroup 
and Markov transition probabilities of Y{t) related by: 


Ptf{x)=[ f{y)nt,x{dy) = E^f{Y{t)) t>0,xe 

JR^ 


/e 


')■ 


Here E^, denotes expectation conditional on y(0) = x G K". Define the domain of 
L to be the set: 


Dom(L) = {f G 


‘) I lim ^(P(/-/) exists in 
i-)-0+ t 


‘)}- 


The Kolmogorov equation associated to the SDE problem (4.1.5) is given by: 
(4.1.6) 


du, 


— {x, t) = Lu{x, t) for all a; e 


and t>0 


with initial condition: 

u{x, 0) = /(x) for all x G K" and for any / G Dom(L) . 


The solution u(x,t) to (4.1.6) has the following stochastic representation: 


u(x, t) = E^f{Y (t)) = Ptf{x) . 

Under suitable conditions on the drift field (e.g. those stated in Assump¬ 


tion 4.1), the process Y(t) admits a unique invariant probability measure H that 
is absolutely continuous with respect to Lebesgue measure on M". Let i^(x) be the 
stationary density associated to H that satisfies 

{L*v)(x) = 0 for all x G K” 

where L* is the (formal) adjoint of the infinitesimal generator L of the SDE solution. 
Hence, 

n(L/) = f Lf{x)v{x)dx = 0 , for all / G Dom(L) . 

For the associated semigroup this invariance implies: 

n(-Pt/) = [ Ptf{x)u{x)dx = [ f{x)u{x)dx = n(/) 

J IR”’ J 

for any t > 0. 


Next we introduce (space) discrete analogs of (4.1.4), (4.1.5), and (4.1.6). 
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4.1.3. Approximation. To introduce the generator of the approximation, let 
{ciliLi be the standard basis on K" and let Sh = {xi} be an inhnite Cartesian grid 
over K” with edge length /i, which we refer to as the spatial step size. The precise 
location of this grid is typically set by the choice of initial condition. Let £°°{Sh) 
(or £°° for short) be the Banach space of all grid functions (p : Sh ^ that are 
bounded in the supremum norm: 


(4.1.7) 


Mloo = sup Iv5(a;)| . 
xeSh 


In the context of the SDE problem (4.1.5), the approximate generator given in (2.5.21 
maps a grid function ip : Sh ^ ^ to another grid function Qp : Sh defined as: 


(4.1.8) 


/ h \ 

Qp{x) = X! ( 2^* ) (t’ ( 2 ^ + ^ei) - p{,x)) 


Here we used the shorthand notation pi = p{x)^ei. We also dropped the 
subscript c appearing in (2.5.2), since the theory in this chapter focuses on proving 


properties of this generator only. For any t > 0, the generator Q induces a Markov 
jump process X{t) on Sh, which can be exactly simulated using the SSA specified 
in Algorithm 

to [naiiM 


2.1[ For basic facts about Markov jump processes, we refer the reader 

mwi\. 


Since our assumptions (given below) permit the drift field entering the reaction 
rates in ( |4.1.8[ ) to be unbounded, Q may not be a bounded operator with respect 
to grid functions in £°°. Nevertheless under a weak dissipativity condition on the 
drift field of (4.1.5) given in Assumption |4.1[ we show that the approximation X{t) 
has a stochastic Lyapunov function. This property is then used to show that the 
local martingale 

t 

(4.1.9) M‘^(t) = p{X{t)) - v7(A(0)) - j Qp{X{s))ds , t > 0 

0 

is a global martingale for all </?: S'/I —>■ M satisfying a growth condition, which is 
developed in §4.3| For functions satisfying this condition, Dynkin’s formula holds: 


C 

(4.1.10) ^xp{X{t)) = p{x) + J KxQp{X{s))ds , t>0. 

0 

This class of functions includes locally Lipschitz continuous functions that have a 
local Lipschitz constant that does not grow faster than a sharp stochastic Lyapunov 
function for the SDF solution. Note that the quadratic variation of is given by 


(4.1.11) 


= / {Qp\X{s)) - 2p{X{s))Qp{X{s))) ds . 


By choosing p{X{t)) in (4.1.9) to be the ith-component of X{t), it can be shown 
that the process X (t) induced by the generator Q satisfies the following stochastic 
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equation: 

(4.1.12) 


dX = fi^{X)dt + V2dM , X{0)=x€Sh 


where we have introduced the drift field with ith component 

for 1 < j < n 

and a local (zero-mean) martingale 

M(t) = (Mi(t),...,M„(t))^eK”. 

The quadratic covariation of the components of A4 satisfies: 


(4.1.13) 




' /■* h 

2 cosh{-^^{X {s)))ds iii=j 
Jo ^ 


0 


otherwise 


for all t > 0. (Since jumps only occur in each coordinate direction, the covariation 
of two distinct components of A4 is zero.) We stress that (4.1.12) is shorthand 


notation for a semimartingale representation of the Markov jump process X(t). 


The noise entering (4.1.12) is purely discontinuous and the drift field is less regular 


than in the original SDE. However, this stochastic equation has the nice feature that 
realizations of X can be produced using the SSA specified in Algorithm |2.1[ This 
stochastic equation resembles an ODE with random impulses 1123]. Properties of 
strong solutions to (4.1.12) are analyzed in (4.3 We also remark that although 
there is a 1/h prefactor in to leading order as expected. 

Note that Q evaluated on the grid Sh can be represented as an infinite matrix 
Q with the following nonzero off-diagonal transition rates from state x G Sh to state 

y G Sh'. 


Q{x, y) = ^ exp ^ ^fi{x)'^{y - x) 


ii\y-x\=h 


or in words, if ?/ is a nearest neighbor of x. The Q-matrix property given in (2.2.1) 
implies that: 

Ql[x,x) = - ^ Qix,y) . 

y&Sh\{x} 

Recall from Algorithm |2.1[ that the diagonal term of Q is the reciprocal of the mean 
holding time of the approximation X(t) at the state x. 

Consider the semi-discrete (discrete in space and continuous in time) Kol¬ 
mogorov equation: 


(4.1.14) 


du^ 

—^{x,t) = Qu^{x,t) for all x G Sh and f > 0 


with initial condition: 


u’^{x,0) = (p(x) for all x G Sh ■ 


Analogous to the continuous case, this solution can be explicitly represented in 
terms of X{t): 


u^{x,t) = E^(p{X{t)) = Ptip{x) 
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where we have introduced P/*, which is the Markov semigroup associated to Q. In 
terms of which, the domain of Q is defined as: 

Dom(Q) = {(/3 S I lim — ip) exists in . 

t^o+ t 

Let denote an invariant probability measure of Q, which we will prove exists 
under certain conditions on the drift field stated below in Assumption |4.1[ In the 
particular case of a gridded state space, 11^ is also unique. In contrast, if the state 
space is gridless, uniqueness is nontrivial to prove because the process may lack 
irreducibility. We address this point in Chapter 

Let v^{x) be the probability density function associated to 11^, which satisfies 
the identity 

{Q*v^){x) = 0 for all x G Sh 
where Q* is the adjoint of Q. This identity implies that 

Il^{Qp) = 0 for all p G Dom((5) . 

Moreover, an invariant measure 11^ satisfies: 

= E = E = n"(^) 

xGSh. xGSh 

for any t > 0. 

Table [T] summarizes this notation. 



process 

generator 

transition 

probability 

semigroup 

invariant 

measure & density 

SDE 

Solution 

Y{t) G K" 

L 

nt,x 

Pt 

n(^) = /a v{x)dx 

A c K” 

Numerical 

Approximation 

X{t) G Sh 

Q 

nt 

TDh 

AcSh 


Table 1. Tabular summary of (space) continuous & discrete objects. 


4.1.4. Drift Field Assumptions. Here we state and illustrate our assump¬ 
tions on the drift field. 


Assumption 4.1. The drift field /i : M” K” of the SDE in ( |4.1.5[ ) is of class 
C"* and satisfies conditions (Al) - (A3). 

(Al) One-Sided Lipschitz Continuity: There exists G K such that: 

Dp{x){^,o<c^\e 

for all X, C e K”. 

(A2) Polynomial Growth: There exist m G No and Cp > 0 such that: 

SUP V SUP 

1 + |xp-+l 1 + |xp-+l-fc - P 

for any 1 < A: < 4. 
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(A3) Weak Dissipativity: There exist /3o > 0 and ao > 0 such that: 

(4.1.15) sign(a:i) < cio (1 + |a:p™) - Po 

for all X S K", for every 1 < i < u, and where m is the natural number 
appearing in Assumption (A2). 

(A4) Minimal Growth: For any £ S N, there exist /3f > 0 and > 0 whose 
growth in £ is subfactorial (i.e. is less than the growth of [£!]) which satisfy: 

> -a. (1 + 

for all X € K”, for every 1 < i < n, and where m is the natural number 
appearing in Assumption (A2). 


We chose these assumptions with the following drift fields in mind. 


Example 1. Let g : K" —?> K" be of class C'^ with component form 

9{x) = 


that satisfies the polynomial growth condition: 


(4.1.16) 


sup 


\9{x)\ ll-P^g(a^)ll 

1 + |x|2™ 1 + 


< Cf, 


for any 1 < A: < 4 


and the one-sided Lipschitz property: 

Dg{x){^, 0 <Cl\e forallx,eGK” 

for some constants Cp > 0 and Cp G M. The drift field whose ith-component can 
be decomposed as: 


fj,i{x) = -Qi xf^^^ + gi{x) , Oi > 0 , 1 < i < 


satisfies Assumption |4.1| Indeed, 

DKx)i^,0= 

l< 2 ,j<n ^ 




where 5ij is the Kronecker delta. This inequality implies that Assumption (Al) 
holds with Cp = Cp. According to the hypotheses made on g, it is apparent that 
Assumption (A2) holds. We also have that 

/i,(x)sign(xi) = -Oj -£ 5 i(x)sign(xj) 

<-a, |x,|2-+i+C®(l + |xp“) 


for all X G M", which implies Assumption (A3) holds with = Cp and /3o = 
mini<i<„ai. By Bernoulli’s inequality and (4.1. there exists ae > 0 such that: 


|M.(x)r = af X 


2e (2m+l)(2e) 


1 - 


9i{x) 

a.x2-+i 


2£ 


> af _ 2^af-^C^p{l + |x|2’”)xf 


for all X G M", which implies Assumption (A4) holds with constants Pi = /3g^ and 
Ui = 2iCp maxi<i<„ whose dependence on £ is subfactorial. 
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Lemma 4.1.1. Assumption\4. 1\ (Al) implies that pL satisfies 


iniy) - fj.{x),y - x) < C^\x-y(^ 


for all x,y € 


Proof. Use Taylor’s formula and Assumption 4.1 (Al) to obtain: 


ih'iy) - Kx),y - x) = / Dy{x + s{y - x)){y - x,y - x)ds 
Jo 


<C^ \y-x\'- 


as required. 


□ 


Remark 4.1.1 (On Usage of Assumption 4.1). Local Lipschitz continuity of 
the drift field is sufficient to show that for any initial condition x G K" and for any 
Brownian motion W there a.s. exists a unique process Y adapted to W that satisfies 


(4.1.5) up to a random explosion time. The weak dissipativity condition (A3) is 


used in Lemma |4.2.3| to prove that Y possesses a stochastic Lyapunov function, 
which implies that this pathwise unique solution can be extended to any t > 0 
[67] . Additionally, Lemma [4. 2.3 is used in Lemma 4.3.1 to obtain bounds on every 
moment of the solution that are uniform in time. 


Assumption (Al) and (A2) are used to prove Lemma 4.3.6 which implies that 
the semigroup Pt is strongly Feller. Local Lipschitz continuity of the drift field is 
also used in Lemma |4.2.4| to prove that the process is irreducible and its transi¬ 
tion probabilities are absolutely continuous with respect to Lebesgue measure. The 
strong Feller property and irreducibility imply uniqueness of an invariant probabil¬ 
ity measure. In fact, a stochastic Lyapunov function, irreducibility, and the strong 
Feller property imply that the process Y is geometrically ergodic with respect to 
this unique invariant probability measure, as stated in Theorem |4.2.5| 

Regarding the stability of the approximation. Assumptions (A2), (A3) and (A4) 
are used in Lemmas |4. 2. 6| and |5. 4.1 to prove that the approximation has a stochastic 
Lyapunov function. In the particular case of gridded state space, the approxima¬ 
tion is also irreducible, and it then follows that the approximation is geometrically 
ergodic with respect to a unique invariant probability measure as stated in Theo- 
Assumption (A2) is used to bound the growth of fi by this stochastic 


4.2.8 


5.2.1 


4.3.11 In the gridless state space context, 
to show that a mollified generator is 


rem 

Lyapunov function in Lemmas |4.3.8 and 
Assumption (A2) is also used in Lemma 
a bounded linear operator in the standard operator norm. This property implies 
that the semigroup of the approximation is Feller, and together with a stochastic 
Lyapunov function, that the approximation has an invariant probability measure. 

Regarding the accuracy of the approximation, we require the additional as¬ 
sumption that the drift field has four derivatives. This assumption is used in 

’■) then 


Lemmas 


4.3.5 


and 


4.3.6 


The latter Lemma implies that if / G 


cti 


Ptf G for any t > 0. This result is then used to estimate the global er¬ 


ror of the approximation over finite (resp. infinite) time intervals in Theorem 4.5.1 
(resp. Theorem 4.5.2). Four derivatives of the drift field are needed for proving 
second-order accuracy because the remainder term in {Qf{x) — Lf{x)) involves 


four derivatives of /(x), see Lemma 4.4.1 


The following Lemma relates our weak dissipativity assumption to more stan¬ 
dard ones. Note that the converse implication does not hold, in general. 
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Lemma 4.1.2. Assumption (AS) implies that there exist constants a > 0 and 
/? > 0 such that 

p{xfx < a (1 + |a;p™+i) - jS |a;p™+2 for all a; € K” . 

Proof. Multiply both sides of ( 4.1.15| ) by |a;i| and sum over i to obtain: 

pi{xfx<ao {l + \x?n\Ai-MAlZXl- 

Since norms on M" are equivalent, there exist Cg, Dq > 0 such that: 

|a^|i < Cq\x\ and ^ Dq\x\'^'^'^‘^ for all x S K" . 

These inequalities combined with the preceding one imply the desired inequality. 

□ 


4.2. Stability by Stochastic Lyapunov Function 


In this section we show that the approximation and the SDE solution are long¬ 
time stable under Assumption |4.1| on the drift field. Specifically, we show that the 


Markov processes induced by the generators L in (4.1.4) and Q in (4.1.8) are geo¬ 


metrically ergodic. The model problem considered for this analysis is (4.1.5), which 


is an SDE of elliptic type with an unbounded state space, locally Lipschitz drift 
coefficients, and additive noise. The weak dissipativity condition Assumption |4.1| 
(A3) is a crucial ingredient to prove that both the SDE solution and numerical 
approximation are geometrically ergodic. We emphasize that for this model prob¬ 
lem the state space of the approximation is gridded and the generator Q can be 
represented by an infinite matrix. Our main tool to prove geometric ergodicity is 
Harris Theorem, which we recall below. We use the total variation (TV) metric to 
quantify the distance between probability measures. 


Definition 4.2.1. If Hi and U 2 are two probability measures on a measurable 
space H, their total variation is defined as: 

(4.2.1) iirii — n2||Tv = 2sup|ni(A) — n2(A)|, 

A 

where the supremum runs over all measurable sets in Q. In particular, the total 
variation distance between two probability measures is two if and only if they are 
mutually singular. 


The factor two in this definition is convenient to include because this metric 
can be equivalently written as 


lini 


n 2 ||TV = sup 



/ if dn 2 s.t. 
Jn 


ip : n —>■ M. , 

\^{x)\ < 1 for all X S D 


4.2.1. Harris Theorem. Harris Theorem states that if a Markov process 
admits a stochastic Lyapunov function such that its sublevel sets are ‘small’, then 
it is geometrically ergodic. More precisely, Harris Theorem applies to any Markov 
process X(t) on a Polish space D, with generator L, Markov semigroup Pt, and 
transition probabilities 

nt_a;(A) = P(x(t) e A I x(o) = x) t>o, xeH, 

that satisfies Assumptions |4.2| and |4.3| given below. 
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Assumption 4.2 (Infinitesimal Drift Condition). There exist a function $ : 
n —)■ K+ and strictly positive constants w and k, such that 

(4.2.2) L$(a:) < k — w <i)(a:) , 

for all a; € D. 


Assumption |4.2| implies that the associated semigroup P* satisfies: 

Pi$(a:) < e-™*$(a;) + -(1 - e"™*) 
w 

for every t > 0 and for every x G M". 

Assumption 4.3 (Associated ‘Small Set’ Condition). There exists a constant 
a S (0,1) such that the sublevel set {z G | $( 2 ;) < 2k/w} is ‘small’ i.e. 

(4.2.3) ||nt,,-nt,J|Tv<2(l-a) 

for all x,y satisfying $(a:) V $(?/) < 2k/w, where the constants k and w are taken 
from Assumption |4.2[ 

Remark 4.2.1. Assumption |4.3| is typically verified by proving that: 

(1) There exists an accessible point: that is a 1 / G such that 

(4.2.4) ^tAOy) > 0 

for all X G n, for every neighborhood Oy of y, and for some t > 0. 

(2) For every t > 0 the associated semigroup is strongly Feller. 

(4.2.5) G Bb{n) G ^(r!) . 


If a Markov process has an accessible point and its semigroup is strongly Feller, 
then every compact set is small (including the sublevel set of the Lyapunov function 
mentioned in Assumption 4.31 and the process can have at most one invariant 
probability measure [I06llim l43]. 


Theorem 4.2.1 (Harris Theorem). Consider a Markov process X(t) on H 
with generator L and transition probabilities P\t,x, which satisfies Assumptions 4^.2 
and 4-3 Then X{t) possesses a unique invariant probability measure PI, and there 
exist positive constants C and r (both depending only on the constants w, k and a 
appearing in the assumptions) such that 

||nt,x - PIIItv < C exp(-r t) $(x) , 


for all X G Cl and for any t > 0. 


Remark 4.2.2. It follows from (4.2.2) that: 


n(4>) < - . 

w 


For further exposition and a proof of Harris Theorem in a general context, see 
the monograph [106j (or |46j for an alternative proof), and in the context of SDEs, 
see Section 2 and Appendix A of [53j . On a countable state space, the conditions 
in Theorem |4. 2. 1| become less onerous, and in fact, if the process is irreducible then 
all finite sets are small, in the sense of Assumption |4.3| In fact, it is sufficient to 
check irreducibility and an infinitesimal drift condition una see Theorem 7.1]. 
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Theorem 4.2.2 (Harris Theorem on a Countable State Space). Consider an 
irreducible Markov process X on a countable state space VL with generator L and 
transition probabilities which satisfies the infinitesimal drift condition As¬ 

sumption f.2. Then X possesses a unique invariant probability measure □, and 
there exist positive constants C and r (both depending only on the constants w and 
k appearing in the assumption) such that 

lirit^x - HIItv < C exp(-r t) $(a;) , 
for all X G n and for any t > 0. 

4.2.2. Geometric Ergodicity of SDE Solution. The material here is an 
extension of Theorem 2.12 in m to non-self-adjoint SDEs. We first show that a 
weak dissipativity condition on the drift implies that the generator L of the SDE 
(4.1.51 satisfies an infinitesimal drift condition. 


Lemma 4.2.3. Suppose the drift field of the SDE (4.1.5) satisfies Assump- 
tion\4-l\ Consider the function V : K" —?► 


(4.2.6) 


V{x) = exp (a|x|2'"+2) 


+ defined as: 
for any 0 < a < 




2(m -I- 1) 


where m and (3 are the constants appearing in Assumption 4.1 (A2) and Lemma 4-1.2. 
respectively. Then there exist positive constants K and 7 such that 

(4.2.7) LV{x) < K-jV(x) 

holds for all x G K". 

Remark 4.2.3. Note that the Lyapunov function in ( |4.2.6 ) is defined in terms 
of the natural number m appearing in Assumption |4.1| (A2), which sets the leading 
order polynomial growth of the drift field fi. This Lyapunov function is sharp in 
the sense that it implies integrability of any observable that can be dominated by 
the exponential of the magnitude of the drift field. A discrete analog of this lemma 
is given in Lemma |4.2.6[ and is a crucial ingredient to the subsequent analysis of 
the approximation. 

Proof. Set £{x) = x^x and 0 (m) = exp(aM*’), so that V = Qo£. We assume 
that p is a free parameter at first. The proof then determines p such that the desired 
infinitesimal drift condition holds. Erom these definitions it follows that: 

D£{x)^r] = 2x'^ri , D'^£{x)(q,ri) = 2r]T]^ , 

0'(m) = apu^~^Q{u) , and 0 "(u) = apfp — 1 )m^“^)0(m) . 

Hence, 

L{Qo£) 1 


Qo£ 


0o£ 

1 


0o£ 

= ap£P-^L£ 


{{O' o £)L£ + {&' o £)\D£\‘^) 

{{& o £)L£ + o £)£) 

Aa‘^p^£‘^P-^+Aap{p-l)£P-^ 


Lemma 4.1.2| implies that 

L£ < -2/3£'^+^ + 2 a£™+i /2 + 2{a + n) . 
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Applying this inequality to the preceding equation yields, 

L(0o£) 


0o£ 


< 2ap(2ap£^P-^ - /3£P+”^) 

+ (4ap(p — 1) + 2ap{a + n) + 2a£™‘^^^^)£P~^ . 


Given the hypothesis on a, the negative term —j3 £p+'^ in this upper bound dom¬ 
inates over the leading order positive term 2 a p £^p~^ provided that p < to -|- 1 . 
Hence, the largest p can be is to -I- 1. Substituting this value of p yields: 


L(0o£:) 

0 o£ 


< 2a{m + 1 ) {2a{m + 1) — fi) -p lower order in £ terms . 


< 0 by hypothesis 

Hence, for any 7 > 0, there exists E* such that: 
L{Qo£) 


eo£ 


< —7 for all \xf > E* . 


Accordingly, a suitable constant K in ( 4.2. 7| ) is: 

K= sup |L(0 o £)-1- 7(0 o £)| 

|x|2<E* 

which is bounded since the set {|a:p < E*} is compact and 0 o £ is smooth. □ 

Lemma 4.2.4. Suppose the drift field of the SDE ( |4.1.5 ) satisfies Assump¬ 
tion^^ For every t > 0 and E > t), there exists e G (0,1) such that 

(4.2.8) ||nt,,-nt,J|Tv<2(l-e), 

for all x,y satisfying V{x) V V{y) < E. 

Proof. It follows from the assumed ellipticity of the equations and the as¬ 
sumed regularity of the drift field, that there exists a function p(t,x,y) continuous 
in all of its arguments (for t > 0) such that the transition probabilities are given 
by nt^x(dy) = p(t, x, y) dy. Furthermore, p is strictly positive (see, e.g.. Prop. 2.4.1 
of |17p. Hence, by the compactness of the set {x : V{x) < E}, one can find a 
probability measure f] on K" and a constant e > 0 such that, 

for any x satisfying V{x) < E. This condition implies the following transition 
probability Pt on is well-defined: 

1 


H*,; - 


-41t X — 




1 — e 1 — e 

for any x satisfying V{x) < E. Therefore, 

— Ht^j^llxv = (1 — e)l!nt^2, — Ht^yllxv 

for all a:, y satisfying V(x) \/ V(y) < E. Since the TV norm is always bounded by 
2 , one obtains the desired result. □ 

With Lemmas |4.2.3| and |4.2.4[ we can now invoke Harris Theorem |4.2.1| to 
obtain the following geometric ergodicity result for the SDE solution. 
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Theorem 4.2.5. Suppose the drift field of the SDE (4.1.5) satisfies Assump¬ 
tion 4-1 Let V{x) be the Lyapunov function from Lemma 4-2.3. Then, the SDE 


admits a unique invariant probability measure 11 . 
constants A and C such that 


Moreover, there exist positive 


(4.2.9) 


\^t,x - n||Tv < cexp(-Xt)V{x) 


for all t > 0 and for all x G K". 

Proof. According to Lemma 4.2.3 for any a G (0,/3(2(m + 1))“^), V(x) is a 
Lyapunov function for the SDE. Moreover, by Lemma [4. 2. 4| it satisfies a small set 
condition on every sublevel set of V. Hence, Harris Theorem implies that the SDE 
solution has a unique invariant probability measure H and that (4.2.9) holds. □ 


4.2.3. Geometric Ergodicity of Approximation. Here we show that the 


Markov process induced by the generator Q in (4.1.8) is geometrically ergodic. The 


main ingredient in this proof is a discrete analog of Lemma |4.2.3| which shows that 
the approximation has a stochastic Lyapunov function. 


(4.2.10) 


Lemma 4.2.6. Suppose the drift field of the SDE satisfies Assumption \4-l\ Let 

Po 


V(x) = exp {a\x\lZXl) for any Q < a < 


2{m + 1 ) 


where m and /3o are the constants appearing in Assumption 4-1 (A2). Then, there 
exist positive constants h^, K and 7 (uniform in h) such that 


(4.2.11) 

holds for all x G K" and h < he 


QV(x) < K — 7 H(x) 


Remark 4.2.4. Note that the stochastic Lyapunov function in (4.2.10) dit 


fers from the one defined in equation (4.2.6) of Theorem 4.2.3 In particular, the 


Lyapunov function in (4.2.6) uses the 2-norm whereas the one in (4.2.10) uses the 
(2m -b 2)-norm. Using the (2m -|- 2)-norm in (4.2.10) is convenient for the analysis 


and makes the proof more transparent, since the reaction channels in Q are the 
coordinate axes. However, since all norms in K" are equivalent, this difference is 
minor, and in particular, does not affect the type of observables we can prove are 
integrable. 


The following lemma is used in the proof of Lemma |4.2.6| 

Lemma 4.2.7. Suppose that a, b > 0 and k G N. Then we have the bound: 


(4.2.12) 


- < k (a - b) a 


k-l 


Proof. Note that the bound holds if a = 0 for any b > 0 and for any k G N. 
For a > 0, the bound comes from Bernoulli’s inequality 

k , 


(4.2.13) 


>l + k(--l). 


Indeed, recall that Bernoulli’s inequality holds for all b/a > 0 and k € N, with 
equality if and only if a = b or k = 1. Since a > 0, (4.2.13) implies (4.2.12). □ 


Here is the proof of Lemma |4.2.6[ 
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Proof. Apply the generator Q in ( |4.1.8[ ) to V{x) to obtain: 
(4.2.14) 

where we have defined: 




for any 1 < i < n. 

Step 1: Bound Af. Lemma 


vtM 

V{x) 


-1) , (x) = V{x ± e,h) 


4.2.7 


implies that: 


with 

a = {xi± h)^ , h = xf , and k = m + 1 


^ = exp ( a((x. ± hfr+^ - a{xjr+^ 


< exp ^o(to + l){±.2xih + h^){xi ± 

< exp f ± 2a(m + l)x^'^+^h + Rf 


where we have introduced: 

Rf = a{m + l)x1'^h'^ + a{m + l)(±2aiih + h?) ({{xi ± — (a;^)'") . 

A second application of Lemma |4.2.7| shows that 

liix, ± hfr - < m\ix. ± h)^ - V ix, ± hfr-^ 

from which it follows that 

Rf < a{m + l)xf^h^ + am{m + l){±2xih + h^Y{2x1 + 

< a{m + l)xf^h^ + am{m + l){8x'fh^ + 2h'^){2xj + 2h?)'^~^ 

< {Cixj^ + Co)h2 

for some CojCi > 0, for all x G K" and h < 1. Thus, we obtain the following 
bound: 


(4.2.15) 


Af < ^ exp ( ±fp,i 


X ^ exp (±2a(m + \)xf^^^h + {Cixf^ + Co)h^) — 1^ 


Step 2: Merge bounds. Substitute (4.2.15[) into (4.2.14) to obtain: 


(4.2.16) 


QVjx) 

Vix) 


n 

< (l){h,x) = y^(/)i(h,a:) 


where we have introduced: 

2 

(4.2.17) (f)i{s,x) = (exp(ci(a:)s^) cosh(ai(a;)s) — cosh(6i(x)s)) 
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and the following functions of position only: 

ic,{x)=Cixj^ + Co 
( 4 - 2 - 18 ) < bi{x) = 

[ai(x) = 2a(m + + bi{x) 


Step 3: Leverage monotonicity. We propose to upper bound the right-hand-side 
of (4.2.16) by the limit of (f){h,x) as h —>■ 0. Indeed, by rHopital’s Rule this limit 
is given by: 


lim x) = ^2 ~ + 2ci 

i=l 

n 

= ^ 2a{m + l)x‘f'^{2a{m + + fJ-iXi) + 2ci 

71 

= ^ 2a{m + l)(2a(m -I- 1) - H- 

i=l 

(4.2.19) < 2a{m + 1) (2a(m -I- 1) — Po) + lower order in |x| terms 

< 0 by hypothesis on a 

where in the third step we used Assumption |4.1| (A3) to bound fiiXi = fXi sign(a:i)|a:i|. 
Thus, this upper bound can be made arbitrarily large and negative if \x\ is large 
enough. Another application of rHopital’s Rule shows that 

(4.2.20) lim^(s,a;) = 0 

s^o ds 

and hence, the point s = 0 is a critical point of 4>{s, x). To determine the properties 
of this critical point, we apply THopital’s Rule and Assumption |4.1| (A4) to obtain: 

1™ = ^—{at-bt + 12a-c, -b 12c-) , 

n 

n 

< ^2 + l)(2a(TO -b 1) — l3o)x^"^^'^{a'f + bf) + ■ ■ ■ 

i—1 

71 ^ 

- 55 + l)(2a(m -b 1) - + ■ ■ ■ 

^ O 

- 55 + l)(2a(m -b 1) - + • • • 

2^1 

(4.2.21) < ^a(m + 1) (2a(m + 1) — /?o) kl 8 m +4 + lower order in \x\ terms 

24 ^ ^ ^ 

< 0 by hypothesis on a 


which is negative if |a;| is large enough. By the point-wise second derivative rule, 
lim^-j-o </>(s,a;) is a local maximum of the function ^(s,x), and therefore, the func¬ 
tion (p{s, x) is strictly decreasing in a neighborhood of s = 0. Thus, we may use 
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lims_>o in (4.2.191 to bound (j){s,x) as long as s is sufficiently small (in order 

to also justify the bound on Af in Step 1) and |a;| is large enough. In the same 
way, we can bound a Maclaurin series expansion of ^(s, x) as follows 

0„2(fc-l) 

0 (s, x) - lim cj,is, x)= {2k)\ + 


l<i<n 

2<fc<oo 


= E 


l<i<n 

2<fc<oo 


^ E 


l<i<n 

2</c<oo 


2s2(fc-i) 

2s2(fc-l) 


(a? - E 


2(fc-,) 20-1) 




2 a(m + l)( 2 a(m + 1 ) - Po)x‘l'^^'^ | ^ j + • • • 

i=i 


^ E 

i<i< 

2<k<< 

^ E 


l< 2 <n 

2<fc<oo 


2(s/2)2('=-i) 

m'- * 

2(s/2)2('=-i) 

M'- 


2 a(TO + l)( 2 a(m + 1 ) -H- 

2/3fc_ia(m + l)(2a(m + 1) - H- 


2(s/2)2('" , 11 /r, ^ ,11 oil |(4m+2)fe 

- E --2/3fc_ia(m + l) {2a{m + l)-0o) 

< 0 by hypothesis on a 

+ lower order in I a: I terms 


Since Assumption 4.1 (A4) states that the constants j3e. and ai are subfactorial in 
G N, the series appearing in this bound are convergent. This bound shows that 
the leading term in |a:| at every order in the Maclaurin series expansion is negative, 
and therefore, the upper bound lims_io (t>{s^ x) on (()(s, x) holds with distance from 
the origin. 

Step 4- Conclude Proof. We have shown that there exist positive constants 
R* and 7 (all independent of h) such that: 


(4.2.22) 


QV{x) 

V{x) 


< —7 for all a; G K" satisfying \x\ > R* 


and for all h < he- Define K in (4.2.11) to be: 

K= sup |LI4(a;) + 7 D(a;)| + sup sup \{Q — L)V[x)\. 

|2,|<i?* 0<h<h^ |a:|<fl:* 

The first term of K is independent of h and bounded, since LV(x) + jV(x) is 
smooth and the set {|a:| < R*} is compact. According to Lemma 4.4.1 below, the 
difference \{Q — L)V{x)\ is 0{hf), involves four derivatives of V{x) and the drift 
field pL. Since this function is continuous and the set {|a:| < R*} is compact, the 
last term can be bounded uniformly in h. □ 

We are now in position to invoke the countable state space version of Harris 
Theorem |4. 2. 2[ which implies geometric ergodicity of the approximation. 
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Theorem 4.2.8 (Geometric Ergodicity of Approximation). Suppose the drift 

Then the Markov 


field of the underlying SDE (4.1.5) satisfies Assumption f.l 


process induced by Q has a unique invariant probability measure supported on 
the grid Sh ■ Furthermore, there exist positive constants A and C such that 


— Xt 


(4.2.23) ||nt-n^||Tv<GE(a:)e 

for all t > 0 and all x € Su- 

Recall that on a countable state space the TV distance between probability 
measures is equivalent to the £i distance between their probability densities. Thus, 


(4.2.23) is a bound on the ii distance between the transition and stationary prob¬ 


ability densities of the approximation. 

Proof. The Markov process induced by Q is irreducible. Indeed, at every 
point on the grid, the transition rate to nearest neighbors is strictly positive. Thus, 
for any two grid points x and y, it is possible to find a finite sequence of jumps 
from X to y, and for any t > 0 , there is a strictly positive probability that this finite 
sequence of jumps will occur. Moreover, the generator Q satisfies an infinitesimal 


drift condition according to Lemma 4.2.6 Hence, Theorem 4.2.2 implies that (i) 


the process possesses a stationary distribution, and (ii) the transition probability 
of the process converges to this stationary distribution geometrically fast in the TV 
metric. □ 

Remark 4.2.5. Let || • ||y denote the V-weighted supremum norm defined as: 
WfWv = sup ^ 

where / : S/j —>■ K and V : K” —>■ K is a Lyapunov function of Q from Lemma [4.2.6| 
Define the following set of grid functions: 

£v = {g ■ Sh I ||5||v<oo}. 

The conditions of Theorem |4.2.8| imply exponential convergence in this function 
space: 

(4.2.24) ||P,V - n\f)\\v < C\\f - n'‘(/)||ve-^‘ , 

for some A > 0 and C > 0, for alH > 0 and for all f G £v [TTsl [731174]. Here we 
used the notation: 

n"(.f) = E ■ 

xeSh 

Remark 4.2.6. Consider the Hilbert space 

= {f:Sh^R\ Y. 

xeSh 

with inner product: 

(/,ff)i;2(nD = Y fi^)gi^)'^^i^) ■ 

xGS 

If Q is self-adjoint, in the sense that 

if, Qg)e^{n>-) = {Qf, 9)i^n>') for all f,g G 4(n'*) 
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then (4.2.23) holds if and only if there exists a A > 0 such that: 
(4.2.25) |1P,V - n^/)||, 2 (n.) < C\\f - n^/)||,.(n.)e-^‘ 


for alH > 0 and for all / G t' 2 (n^) [1181118) . We stress that this equivalence does 
not hold in general if the Markov process is not self-adjoint. 


Remark 4.2.7. Let (j{Q) denote the spectrum of Q: the set of all points A G C 
for which Q — A / is not invertible. If Q possesses a spectral gap in i.e., a 

gap between its zero eigenvalue and the eigenvalue of second smallest magnitude, 
then its transition probabilities converge exponentially fast in £^(n^). Indeed, by 
the spectral mapping theorem, 

= I AGa(g)} 

and hence, a spectral gap in Q implies that has a gap between its unit eigenvalue 
and its second largest eigenvalue [mill]. 


4.3. Properties of Realizations 


4.3.1. Basic Properties of SDE Solution. Here we review some basic prop¬ 
erties of the SDE solution under Assumption |4.1| on the drift field. 


Lemma 4.3.1. Suppose the drift field of the SDE (4.1.5) satisfies Assump¬ 
tion \4 A\ For every ip : M" - 
(/? € Dom(L), we have that: 


satisfying \ip\ < V + Cq for some Cq > 0 and 


\EMym<v{x) + - + co 

7 

for all t > 0 and for all x G K". Here K and 7 are taken from Lemma 4-2.3\ 

As a consequence of this lemma, all moments of the SDE solution are bounded 
uniformly in time. 

Proof. Since (p is dominated by V, this result follows directly from Lemma [4.2.3| 
Indeed, 

K 


\Pt^ix)\ < Pt\v>{x)\ < PtV{x) + Co<e yV{x) -\ -1- Co . 

7 


□ 


The next lemma bounds the mean-squared displacement of the SDE solution. 


Lemma 4.3.2. Suppose the drift field of the SDE (4.1.5) satisfies Assump- 
tion \4-l\ Then the displacement of the SDE solution Y satisfies: 

Ex\Y(t) — x\'^ < 2 t"^ ^(x)C q^ 4: n t 

for all t >0 and for all x G M". Here K, 7 and Cq are taken from Lemma \4.3.1\ 
Proof. Apply Jensen and Cauchy-Schwarz inequalities to obtain 
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Take expectations, use the polynomial growth of fj,, apply Lemma |4.3.1[ and invoke 
the Ito isometry to derive the desired result. □ 


The one-sided Lipschitz property implies that the SDE solution is almost surely 
Lipschitz with respect to its initial condition as the next lemma states. 


Lemma 4.3.3. Suppose the drift field of the SDE (4.1.5) satisfies Assump¬ 


tion 4-1 For s < t, let Yt s{x) denote the evolution map of the SDE, which satisfies 
Ys^s(x) = X and Yt^s o ys,r{x) = Yt^r(x) for r < s < t. This evolution map satisfies: 

\Ys+t,s{x) - Ys+t,s{y)\'^ < |a; - y\^ exp(2 C7^ t) a.s. 

for all t > s and for all x,y G M". Here is the one-sided Lipschitz constant of 
the drift field from Assumption 4-1 (-^V- 

Proof. By Ito-Taylor’s formula: 

\Ys+t,six) - Ys+t,siy)\'^ 

— 1^ y\ y I (^-t-r,s(^) (y))/^(Ls+r.s (^)) (?/))) 

Jo 

<\x-y\‘^ + 2Cf f \Ys+r,six) - Ys+r,siy)\‘^ dr 

Jo 

where we used the one-sided Lipschitz property of the drift field with Lipschitz 
constant Cf. Gronwall’s Lemma implies the desired result. □ 

We now turn to estimates related to the semigroup of the SDE solution. 

Lemma 4.3.4. Suppose Assumption \4-l\ holds. Then there exists c > 0 such 

that 


for all x,yG 


\PMx) - Ptipiy)\ < c (1 A t) ||(/j||oo \y - x\ 
^, y} G and t > 0. 


We stress that the test function in this Lemma is not necessarily continuous. 
Assuming that t > 0, this Lemma shows that P* regularizes test functions, and 
specifically that the function Pt(p is Lipschitz continuous for any (p G BbiW^). 
Since Pt<p is also bounded, it follows that the semigroup Pt is strongly Feller, see 
Remark |4.2.1| for a definition of the strong Feller property. This property of P* 
relies on the generator L of the SDE being of elliptic type. For a detailed proof 
of this result that only requires local Lipschitz continuity and a polynomial growth 
condition on the drift field, we refer to §2.3.1 of m- 

Now we make use of the fact that the drift field is four times differentiable and 
these derivatives satisfy a polynomial growth condition. 


Lemma 4.3.5. Suppose Assumption \4-l\ holds. For any p > 1, for any 1 < 
j < 4, for any t > s, and for all x G M", the evolution map Yt^six) is four times 
mean-square differentiable and 

sup E\\D^Yt^six)\\^ < Cpit) IIpIIoo 
where Cp(t) > Q is a continuous increasing function. 


See Theorem 1.3.6 of |17j for a detailed proof of Lemma 4.3.5 under the stated 
assumptions. We apply this result in the next lemma. 
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Lemma 4.3.6. Suppose Assumption \4-l\ holds. For any 1 < j < 4, for any 
t > 0, and for any (p G C'^(R"), the function Pt<p is j times differentiable and there 
exists C{t) >0 such that 


where 


WPtFh < C{t) 
is the norm introduced in (4.1.3). 


In short, this lemma states that Pt maps (M") into itself for any t > 0. 


Proof. For any 1 < j < 4, differentiate j times with respect to x the proba¬ 
bilistic representation of the semigroup: 


Ptp{x) = E(/j(r4,o(a^)) 


and use Lemma 4.3.5 to bound these derivatives of Ytp{x). For example, for j = 1 
the chain rule implies that: 


DPt(p{x){ri) = RD(p{Ytfl{x)){DYtfl{x){rii)) 
for any rj G M”, for all x G M” and for all t > 0. Hence, 


lF)Pt(p(x)(i])l < IE£lip(Yt^o(x))(F>Yt^o(x)(r/))l 

< ElPip(Yt^o(x))(DYt^o(x)(?7))l 

<E\\Dy:>{Yt^o{x))\\ ||Drt,o(a;)|| |? 7 | 

< C'(t)|77|(sup \\Dip{x)\\) 


for some C{t) > 0. Here we used Lemma 4.3.5 and the fact that (p G 


this upper bound is uniform in x, we obtain the desired estimate. 


’■). Since 
□ 


Remark 4.3.1. As an alternative to this probabilistic proof of Lemma [4.3.6| 
which is due to mi, an analytic approach to bounding the derivatives of Ptp is 
available in [90117]. This approach is based on approximating Ptp by the solutions 
to Cauchy problems on compact spaces and assumes conditions which are similar 
to Assumption |4.1 1 


4.3.2. A Martingale Problem. The process X induced by Q has the fol¬ 
lowing classification. 

• A is right-continuous in time with left limits (cadlag). 

• A is adapted to the filtration associated to the driving noise. 

• A is a pure jump process, since the process only changes its state by 
jumps. 

• A has the Markov property, since conditional on A(to) = x, the time to 
the next jump is exponentially distributed with parameter that depends 
only on the state x; and, the jump itself has a distribution that only 
depends on the state x too. 

As such the process admits the following semimartingale representation: 


local martingale 


(4.3.1) 


A(t)=A(0)-f PiX{s))ds + M{t) 


The stochastic equation (4.1.12) is shorthand for this integral form of A. Let us 
prove that the local martingale in this representation is a global martingale. The 
next Lemma is useful for this purpose. 
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Lemma 4.3.7. Suppose Assumption 4-1 
stochastic Lyapunov function from Lemma 
the local martingale: 


holds on the drift field, let V(x) be the 
\4-2.(^ and let ip : ^ R. Consider 


t 

M‘^{t) = ip{X{t)) - ip{x) - j Qip{X{s))ds, t > 0 , 

0 


X( 0 ) = a; . 


If if satisfies the growth condition: 

(4.3.2) Qip^{x) — 2ip{x)Qip{x) < V{x) + C for all x S M" 


for some C > 0, then there exists a positive constant he such that M'^(t) is a global 
martingale for all h < be¬ 


lt is not clear how to relax the martingale condition (4.3.21, since it is based 
on the integrability of a sharp stochastic Lyapunov function for the process. 


Proof. The idea of the proof is to show that the local martingale M^{t) 
is dominated by an integrable random variable, namely the maximum process 
M*{t) = supo<s<t |M'^(s)|. (Recall that a local martingale dominated by an inte¬ 
grable random variable is a uniformly integrable global martingale ms main].) 
To show that the maximum process is integrable, we apply the Burkholder-Davis- 
Gundy inequality with parameter p = 1, which implies that there exists a positive 
constant Ci such that 

E,M*(t) < CiE^^{Mv,Mv){t) < Cl^/E,{Mv,Mv){t) 


where we used Jensen’s inequality with respect to the (concave) square root func¬ 


tion. Since the quadratic variation of M^{t) is given by (4.1.11), the hypothesis 


(4.3.2) implies that there exist positive constants C 2 and C 3 such that 


E,M*it) < (P(X(s)) + C 2 )ds < (P(x) + Cs)t 


Here we used Lemma 4.2.6 which implies that V{x) is a stochastic Lyapunov 


function for X. Hence, M*{t) is integrable and M^(t) is a uniformly integrable 
global martingale, for all t > 0. □ 


Thus Dynkin’s formula (4.1.10) holds for any function satisfying the growth 
condition (4.3.2). The next Lemma illustrates that this class of functions includes 
locally Lipschitz functions. 

Lemma 4.3.8. Suppose Assumption \4-l\ holds on the drift field. Suppose ip : 
K" —>■ M is a locally Lipschitz continuous function 

\(p{y) - <p{x)\ < L,p{x,y)\y - x\ V x, y S M” 
with local Lipschitz constant L,p{x,y) that satisfies the following condition 
L^{x,y) < (7^(1 -f exp(|xp™+^) -f exp(|yp""+^)) V x, y G K” 


for some C,p > 0. Then ip satisfies (4.3.2). 

Since the identity map is globally Lipschitz, this Lemma enables studying so¬ 
lutions to the stochastic equation (4.3.1) at any t>0. 
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Proof. Let G‘^{x) = Q(p^{x) — 2ip{x)Qip{x). Use (4.1.8) to write 


(4.3.3) 


G'^(a;) = U ^ exp(^i/i/2) {ip{x + hei) - ip{x)f 


1 2 
+ exp{-fiih/2) {(f{x - hei) - ^{x)) 


Since ip is locally Lipschitz continuous by hypothesis, 

n 

G'^(.x)<'^C^ (l + exp(|a:p’"+^) + exp(|a; + exp(^i/i/2) 

i=l 

+ C^ (1 + exp(|a;p™+^) + exp(|a; - exp{-pih/2) 

71 

<2 Cl ^ (l + exp(|x|^™'''^)) cosh(/ii/i/2) 

|2™+1 )/i/2 ) 


i=l 


< 2 n Cl (1 + expdxp’"'*'^)) cosh(Gp(l + \x\ 

< V(x) + C 2 

where Ci and C 2 are positive constants. Thus, p{x) satisfies the growth condition 
(4.3.21, as required. Note that in the input to the cosh function, we used the bound: 

\hiix)\ < \p{x)\ < Gp(l + |a;p’"+^) V a; S K" 


which follows from Assumption 4.1 (A2). Also in the last step, to obtain the bound 
by the stochastic Lyapunov function introduced in (4.2.10) we used the following 
inequality: 


,|2m+l 


< \x 


2m+2 


< C,\x\iZV2 


which is valid since all norms are equivalent on 


□ 


As a tool to analyze realizations, we quote I’Hopital’s monotonicity principle 
as a lemma (without proof). 

Lemma 4.3.9. Let f,g : [a, &] —>■ ffi fee continuous functions that are differen¬ 
tiable on (a, fe), and suppose that g'{s) 0 for every s G {a,b). If f {s)/g'{s) is 

increasing/decreasing on {a,b), then the function 

. _ f{s) - f{a) 

® g{s) - g{a) 

is increasing/decreasing on {a,b). 

Next we state an estimate on the mean-squared displacement of the process X, 
which is analogous to Lemma [4. 3. 1[ 

Lemma 4.3.10. Suppose Assumption\4-.l\ holds. Then there exist h^. > 0 and 


Cl > 0 such that the solution of the stochastic equation (4.3.1) satisfies: 

EZX{t) - xp < {Cl + V{x)) t e* , 
for any t > 0 and for all h < he- 

As a consequence, if t ~ G{e) then the mean-squared displacement of the 
approximation is 0(e). We emphasize that this lemma is uniform in the spatial 
step size. This uniformity in h follows from the fact that the stochastic Lyapunov 
function from Lemma 4.2.6 satisfies (4.2.11) with constants K and 7 that do not 
depend on h. 
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Proof. Define ip{y) = {y — x)^{y — x). Since 

\Dy}{y) ■ 77 I < 2{\y\ + |x|)|77| for every x,y,r] G K” 

it follows that ip is locally Lipschitz continuous with a local Lipschitz constant that 
satisfies the sufficient conditions in Lemma 4.3.8 Thus, we may conclude that 
Dynkin’s formula holds for p to obtain the bound: 


¥.xp{X{t)) = p{x)f Qp{X{s)) ds 


= E. 


pt pt n 

: / 2{X{s) — x)'^y^{X{s)) ds+ / 'y^2cosh{y,i{X{s))h/2) ds 

Jo Jo 


piX{s))ds + E, |M^X(s))pds 


ft " 


■ E.t- 


^ 2 cosh(/ri(X(s))/i/2) ds 


Here we used the el ement ary inequality 2x"^y < |a;p + |yp. To bound the term 
apply Lemma 4.3.9 with /(s) = (2sinh(s^i/2))^ and g{s) = s^ to obtain 

f'(s) 2y,iCOsh(sy,i/2)smh(syi/2) 

g'{s) s 


Apply Lemma 4.3.9 once more to obtain 


ns) 

9"{s) 


= (cos(s^i/2)^ + sinh(s/rj/2)^) 


which is monotonically increasing with respect to s. Hence, there exists an > 0 
such that: 

< ^(sinh(/i*y))^ 

for all h < he and for any 1 < i < n. This inequality enables us to bound 
uniformly with respect to h. In particular, there exist positive constants Cq and 
Cl such that 

E^p(X{t))< [ E^p{X{s))ds+ [ {Co+E^V{X{s)))ds 


< 


Ea;(p(X(s)) ds + (Cl + V{x)) t 


Here we used the polynomial growth condition on the drift field from Assump¬ 
tion 4.1 (A2) and the stochastic Lyapunov function of the process from Lemma 4.2.6 


Apply Gronwall’s Lemma to this last inequality to complete the proof. 


□ 


Remark 4.3.2. In (4.5 
with respect to Ptp on both 


we show that for any p G C^(M”), 'd^\sh accurate 

finite and infinite time intervals. Since functions in this 
space are Lipschitz continuous. Lemma 4.3.8 implies that any p G C^(K") satisfies 
(4.3.21, and therefore, Dynkin’s formula holds for any p G C^(K"). 


The next lemma quantifies the complexity of the approximation. 
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Lemma 4.3.11. Suppose Assumption 4-1 holds on the drift field. Let N{t) 
denote the average number of jumps of the approximation X (t) on the interval 
[0,i]. Then there exist positive constants he, Ci and C 2 such that 

N{t)<^(y{x) + C2) 

for all t > 0 and h < he- 

This lemma is not a surprising resultj _ it states that on average the number 


of steps of the SSA method in Algorithm 2.1 is about lt/St\, where t is the time 


span of simulation and St is the mean holding time, which is proportional to h^. 
Since each step requires a function evaluation, this lemma provides a crude estimate 


for the total cost of the approximation. Figures [3.8| and 3.30 verify this estimate 
in the context of a log-normal process that uses adaptive mesh refinement and a 
39-dimensional Brownian Dynamics simulation, respectively. 

Proof. Let {tz}^0 ^ sequence of jump times for the process X. The time 

lag between jumps ti+i — U is exponentially distributed with parameter A(a;): 


X{x) = — cosh{p.i{x)h/2) . 


2 = 1 


On each time interval the process is constant, which implies that 

nt c>o 

/ X(X{s))ds = W X{X{ti)){t A U+i -t Xti) . 

•^0 i=o 

Since the drift field p satisfies a polynomial growth condition from Assumpti on |4.1| 
(A2), it is dominated by the stochastic Lyapunov function V{x) from Lemma 4.2.6 
and hence, 

(4.3.4) E, r X{X{s))ds < 4 r ^.V{X(s))ds < ^ {V[x) + -) . 

Jo ^ Jo ^ 7 

Thus, we can invoke Fubini’s theorem to interchange expectation (conditional on 
A(0) = x) and summation to get: 


E, 


pt 00 

/ X{X{s))ds = '^Ea:X{X{ti)){t A U +1 -tf\ti), 

•^0 i =0 


Since the random variable (t A t^+i — t Ati) is independent of X{X{ti)), 


E, 


pt 00 p 

/ X{X{s))ds = (Ea;{X{X{ti)){t Ati+i - t Ati) \ -TrJ 

•to V 


i=0 

00 


= ^E,(A(A(t,))A(A(t,))-'xq<t) 


2=0 

= Nit) . 

Here x is an indicator function, Jy is the natural filtration of the process X{t), and 
we used the law of iterated expectation to write the total expectation up to time 
ti+i in terms of the conditional expectation up to time t = ti using Ft- The desired 
estimate is obtained by combining the last identity with (4.3.4). □ 
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4.4. Generator Accuracy 

The following lemma improves upon Prop. |2.5.1| by giving a remainder term 
estimate. 


Lemma 4.4.1. Suppose Assumption J^.l holds. Let f G and define 

Rf(x) to satisfy 

Rf{x) = Lfix) - Qf{x) . 

Then there exists Cf(x) > 0 such that: 

\Rfix)\ < hfiCf{x) for every a; S M" . 


Proof. Retracing (and strengthening) the proof of Prop. 2.5.1 using Taylor’s 
integral formula (in place of Taylor expansions) yields: 


(4.4.1) 


Qf (^) + + Si + Ci + Df + ) , 


= -Rf{x) 


where we have introduced: 

l4 df 


A, = cosh (l-s)"ds 


8 dx. 


i JO 


Pih 


= cosh s— (l-s)ds 


4 dxf 

d^f 


i Jo 


cosh 


fiih 


ds 


f-i _ ^ 

* “ 6 axf 7o 

= ^exp ^ ^^{x ± scihfil - sf ds 

Applying the elementary inequality: 

cosh(sr(;) < cosh(u;) for any s G [0,1] and lo G K 
and simplifying yields the desired estimate with 


Cf{x) = Cj2 cosh{fj,ih/2) 


(4.4.2) 


X 


df 


dx,. 




d^f 


dx'f 


\Ti\ 


d^f 


dx^ 


sup 

lsl<l 


d^f 


dxf 


(x + sCih) 


where C is a positive constant. 


□ 


4.5. Global Error Analysis 

In this section it is shown that the Markov process induced by the generator Q 


in (4.1.8) is weakly accurate on finite time intervals. Moreover, it is shown that the 
stationary probability density of Q is nearby the true one. These statements are 
made in the context of the model SDE (4.1.5) under Assumption |4.1[ The main tools 
used in this analysis are properties of the SDE solution and the following variation 
of constants formula for a mild solution of an inhomogeneous, linear differential 
equation on a Banach space. 
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Definition 4.5.1. Let A be the generator of a semigroup on a Banach space 
X. For a given initial condition x & X and a function / : K_|_ —>■ X, consider the 
differential equation: 

\v{t) = Av(t) +f{t) fort>0, 


Then the variation of constants formula 

(4.5.2) v{t) = exp{At)x + f exp{A{t — s))f{s)ds 

Jo 


is called the mild solution of (4.5.1). 


For more background on differential equations on Banach spaces see mi In 
the remainder of this chapter, we use (4.5.2) to derive a formula for the global error 


of the approximation. This formula is then estimated to determine the weak and 
long-time accuracy of the approximation. 


4.5.1. Accuracy Theorems. We are now in position to estimate the global 
error of the numerical method. 


Theorem 4.5.1 (Finite-Time Accuracy). Suppose Assumption 4-1 holds. For 
every ip G C^(IR"), there exists a positive constant such that 

(4.5.3) |P,V(^) - PMx)\ < h^ C^,t Vix) 

for any t > 0 and for all x G Sh- 

The main issue in the subsequent proof is that the generators of the SDE 
solution and numerical approximation may be unbounded. 


Proof. Define the global error of the approximation as: 

(4.5.4) e{x, t) = Ptp{x) — Plf(p{x) for any x G St and t > 0 . 

Note that this global error e(x, t) is a mild solution of the following inhomogeneous 
linear differential equation on the Banach space 

(4.5.5) i{x, t) — Q€{x, t) = {L — Q)Pt(p{x) for any x G Sh and t>0 
with initial condition: 

e(x, 0) = 0 for all x G Sh ■ 

From Dehnition |4.5.1[ a mild solution to this equation is given by the variation of 
constants formula: 

(4.5.6) e{x,t) = [ Pt-s{P - Q)Ps‘p{x)ds . 

Jo 

In order to estimate the global error, it helps to express this formula in terms 
of transition probabilities: 
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Apply Lemma [4.4.1| to obtain 

= I / - Q)Ps^iz)ds 

Jo 


< 


° zeSh 

Ilts,xmL-Q)Ps^{z)\ds 


z^Sh 

rt 


< 


f ^t-sAz)Cp,U^)ds 

^GSh 


where Cp^^{x) is given in (4.4.21 with / = Ps<^. Note that there exists C'o(t) > 0 
such that, 

Cp,^{x) < y(a;) + Co{t) 

for all X € M" and for all s G [0,t]. Indeed, from Lemma 


4.3.6 


if (p G C^(R") then 

Psip G C^(R’^) for any s > 0, and by Assumption 4.1 (A2) on the drift field the 
terms that involve the drift field in (4.4.2) are dominated by the Lyapunov function 
V(x). Hence, we obtain 

|e(x,t)| <h^ f + Co(t))ds 


z^Sh 


< t {V{x) + Co{t)) 


for some C'o(t) > 0. 


□ 


In addition, the Markov process generated by Q accurately represents the sta¬ 
tionary distribution of the SDE. 


Theorem 4.5.2 (Stationary Distribution Accuracy). Suppose Assumption 4-1 


holds. For every ip G we have that 


(4.5.7) 

where 




Cip = sup 


0<h<hc 


Cp^ipds 


and Cp^^{x) is given in (4.4.21 with f = Ps<p. 


Proof. Note that this result is not a corollary of Theorem 4.5.1 since the 
latter gives an upper bound on the global error that is increasing with t. In order 
to obtain the desired estimate, we use a slightly different approach that begins with 
averaging the global error with respect to th e invar iant probability measure of the 
approximation H^ (which exists by Theorem 4.2.8): 


zeSh 


|n"(e(t))| = | ^ n\z)eit,z)\ 

[ (L- Q)Ps<pds 

\Jo 


= |n' 

< 

< 


\{L - Q)Ps<p\ds 

rt 


Cp^pds 


'0 
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where Cp^^(x) is given in (4.4.21 with / = Pgip. Here we used: (i) the Fubini 
Theorem to i nterch ange space and time integrals; (ii) the remainder term estimate 
and (hi) the fact that = n^((/?). Passing to the limit 

□ 


4.4.1 


from Lemma 

as t > oo, and using Theorem 4.2.5 gives the desired estimate. 


Remark 4.5.1. One may be able to prove Theorem 4.5.2 by adapting to this 
setting the version of Stein’s method developed in §6.2 of 


Remark 4.5.2. In the context of smooth initial conditions and a smooth drift 
field with at most polynomial growth at infinity, an estimate of the type: 

\\D^Ptip{x)\\<C{l + \x\‘^)e-'^\ yt>0, Vj>l, VxeK”, 

is available in Theorem 3.1 of [134]. Roughly, the proof of this Theorem shows expo¬ 
nential convergence of the first n derivatives of the SDE semigroup Pt in and 

then uses the Sobolev Embedding Theorem to transfer these estimates to 
i.e., trade regularity for integrability. This procedure leads to point wise estimates 
on the derivatives of the semigroup. One may be able to derive similar estimates 
in this context by assuming that the first n derivatives of the drift field exist, but 
we wish to avoid making this assumption. 











CHAPTER 5 


Analysis on Gridless State Spaces 


To study issues that may arise when the noise is multiplicative (or state- 


dependent) or if the spatial step size is adaptive, consider the generator in (4.1.8) 
but, with random spatial step size equal to where ^ G 17(1/2,1). Figure 5.1 
illustrates the difference between points spaced in this variable way and equally 
spaced points on a grid. In particular, with random spacing the state space of the 
approximation becomes gridless. In this chapter we adapt the stability/accuracy 
results from Chapter to this setting. See §2.2| for the definition of a gridless state 
space. 



(a) grid 



(b) gridless 


Figure 5.1. Gridded vs Gridless State Spaces. This figure illus¬ 
trates the difference between gridded and gridless state spaces. Points 
in these spaces are produced from a given seed in R using fixed spatial 
step size h (left panel) and random spatial step size (right panel) 
where ^ ~ 17(1/2,1). The set of all such states is the state space of the 
Markov jump process. In the case of a gridded state space, every pair 
of points is connected by a finite path on the grid. Thus, if the jump 
rates are positive at every point, then the Markov process is irreducible 
on this gridded state space. In contrast, in the gridless case, there exist 
pairs of nodes such that no finite path on the graph has those nodes as 
endpoints. Thus, even if the jump rates are positive at every point, the 
Markov process is not irreducible on this collection of points. 


5.1. A Random Walk in a Random Environment 

Here we introduce a realizable discretization that uses variable step sizes. To 
introduce this discretization, in addition to the notation presented in Chapter 
let I : M” —)• [1/2,1] with f(x) ~ 17(1/2,1) for all x G K”, i.e., f assigns to each 
point X G K" a uniformly distributed random variable on the interval [1/2,1]. To 
ensure that the coefficients of the resulting generator are continuous, we regularize 
this function by convolving it with the standard mollifier p on K” to obtain: 

(5.1.1) ^{x) = f p{x — y)i{y)dy for all a; G K" . 

/r" 
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Note that 

^min ^ ^{x)h < /imax for all X S K" 

where /imin = h/2 and /imax = h. The statements in this chapter are conditional 
on In particular, we do not average over the random function 

Let r G No- Below we derive conditions on this parameter in order to guaran¬ 
tee that the space-discrete operator is bounded and that the approximation has a 
stochastic Lyapunov function. Given the field the spatial step size parameter h, 
and a test function / : M" —>■ K, define the second-order generator: 



This generator induces an approximation on a gridless state space, which is a 
continuous-time random walk in a random environment |TMII145| . Since a sum 
of discrete random variables is discrete, the state space of this approximation is 
countable. However, the process is not irreducible on this gridless state space as 
illustrated in Figure |5.1[ In addition, if we embed the state space of the approx¬ 
imation into M", the process is still not irreducible with respect to the standard 
topology on K", since at any finite time there may be zero probability for the pro¬ 
cess to reach an open neighborhood of a designated point in state space starting 
from any other point. Even though the approximation may lack irreducibility in 
the standard sense, it turns out the approximation is still long-time stable, weakly 
accurate, and accurate with respect to equilibrium expectations. 


5.2. Feller Property 


Here we prove a regularity property of the semigroup associated to the gridless 
generator Qq. 


Lemma 5.2.1. Suppose Assumption 4-1 holds. If the exponent r in (5.1.2) 
satisfies r > m + 1, then for any h > 0 the linear operator Qq is bounded. 

Proof. It suffices to show that there exists C > 0 such that 


(5.2.1) llQ,/||oo<G 

for an arbitrary f G Bh{ 
reaction rate from x to the state yf{x) = x ± ei^{x)h. By Assumption 4.1 (A2) 


^). In order to prove this statement, let jfj x) b e the 


the drift field satisfies a polynomial growth condition with parameter Cp, and since 
r > TO -f 1, we have that 


< C{h) 
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for all X € K" and for any 1 < i < n. Apply this bound on the reaction rates to 
obtain the following bound on Qqf, 

n 

- fi^)) + - f{x)) 

n 

< 2||/||oo sup {x) + J- (x)) 

<AnC{h) ll/lloo 

which shows that the generator Qq is bounded. □ 

Since Qq is a bounded linear operator, the semigroup it induces is Feller, as 
the next Lemma states. 

Lemma 5.2.2. Suppose Assumption \4-l\ holds. For any t > 0 and r > m + 1 , 
the semigroup generated by Qq is Feller, which means that for each / G 
we have that Pjff G Ct,(R"). 

Proof. Since Qq is a bounded linear operator on the Banach space 
then the operator Pf = exp{Qqt) is bounded on Cb{R'^), and the homogeneous 
Cauchy problem associated to Qq has a unique classical solution for every initial 
condition in Finally, this solution belongs to Cb{RF) for all t > 0. □ 


IIQ9/II00 = sup 

a:GR’* 


5.3. Generator Accuracy 


Here we state that Qq is accurate with respect to the infinitesimal generator of 
the SDE. 


Lemma 5.3.1. Suppose Assumption 4-1 holds. Let / G and set 

Rqf{x) = Lf{x) - Qqf{x) . 

Then there exists C / (x) > 0 uniform in f) such that: 

\Rqfix)\ < hf Cf(x) for every x € M" . 

Proof. This result follows immediately from Lemma [4.4.1 since applying the 
bound from that lemma with h replaced by fh, and bounding the mollifier in Qqf{x) 
by unity, yields: 

\Rqf{x)\<eh^Cf{x)<hl,^Cf{x) 

since fh < h^i,^- Note that the bounding positive function Cf{x) is the same as 
the one appearing in Lemma 4.4. 1| 


Cf{x) = CY^ cosh{iJ,ih/2) 


(5.3.1) 


2=1 


X \Fz 


df 


dxi 




d^f 


dxj 


I Mi 


53 / 


9x? 


sup 

|s|<l 




dxj 


(x + scih) 


where C is a positive constant. 


□ 
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5.4. Stability by Stochastic Lyapunov Function 


Here we show that the approximation on a gridless state space has a stochastic 
Lyapunov function. 


Lemma 5.4.1. Suppose the drift field of the SDE (4.1.5) satisfies Assump- 
tion\4.1\ Let 


(5.4.1) 


V{x) = exp (a|a:| 2 ™+ 2 ) 0 < a < 


Po 


2{m + 1 ) 


where m and Pq are the constants appearing in Assumption 4.1 (A2). Then, there 
exist positive constants K and 7 (uniform in h and f) such that 

(5.4.2) QqV{x) < K--fV(x) 

holds for all x € K", h sufficiently small, and r < 2m + 1 . 

Proof. The argument in this proof is the similar to the gridded state space 
context, see Lemma [4. 2. 6 [ Specifically, the proof shows that there exist he, R* and 
7 > 0 such that: 

QqV{x) 


(5.4.3) 


V{x) 


< —7 for all a: S K" satisfying |a:| > R* 


and for all h < he- We then choose K in (5.4.2) to satisfy: 

K> sup \LV{x)-\-{x)\ + sup sup \{Qq — L)V{x)\ . 

The first two terms of K are independent of h and ^ (x) since V(x) is smooth 
and the set {|a:| < R*} is compact. According to Lemma 5.3. 1[ the difference 
\{Qq — L)V{x)\ is 0(/imax)) iuvolves four derivatives of V{x), and is independent of 
the random environment. Since V{x) is smooth and the set {|x| < i?*} is compact, 
the last term can be bounded uniformly in h. Thus, we obtain the desired estimate 
if we can find R* (uniform in h and f) such that (5.4.3) holds. 

Apply the generator Qq in (|5.1.2[) to V (x) to obtain: 


(5.4.4) 


QqVjx) 

Vix) 


^ Ai(x)exp (-^(x)2 h"^ Ixp’') 


2=1 


where we have defined: 


Ai{x) = 


1 


1 


exp 1 y/li 




^h 

exp I —TrSi 


exp |^a((xi + 5 / 1 )^)""+^ - o((Xj)^)’”+^ 
exp (a{{x^ - ^h)^)'"+^ - a((xi)^)'"+^ 


- 1 


- 1 


Note that Ai involves positive and negative terms, and in order to obtain the 
estimate (5.4.3) we must show that the negative terms dominate in a way that is 
uniform in h and f. For any x G K", note that Lemma[4.2.7|implies that Ai satisfies 


. / N 1 


1 


for * = 1 ,..., n. 


exp 


fh 


fh 




exp \^{2xifh + + l)(xi + 

exp (a{—2xifh + + l)(xi — 


- 1 


- 1 
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These bounds on Ai imply that: 
(5.4.5) 


< (/'(/i, a:) = ^ (pi(h, x) 

V(x) 


where we have introduced: 

(5.4.6) 

2 

4 > i { s , x ) = (exp(cj(ai)^^s^)cosh(ai(a;)^s) - cosh( 6 i(x)^s)) exp {-^‘^s‘^d{x)) 
and the following functions of position only: 


(5.4.7) 


ai{x) = 2 a (m + l) 

h{x) = 

c,(x) =Cixf™ + Co 
d{x) = \x 


\2r 


where Ci and Cq are positive constants. 

We now derive an upper bound on the right-hand-side of (5.4.51 from the limit 
of 4'{h,x) as h —>• 0. By rHopital’s Rule this limit is given by: 

r j.f \ (®i ~ 

hm^(s,x)=> " I' 

i=l 

(5.4.8) = 2a(m -I- l)(2a(m -f 1) — P)\x\'4mXi + lo'^^r order in \x\ terms 

Note that if |x| is large enough, lims_j.o (('(s, a:) is large and negative. Another 
application of I’Hopital’s Rule shows that 

d(j) 


(5.4.9) 


lim T^(s, x) = 0 
s-s-O os 


and hence, the point s = 0 is a critical point of ^(s, x). To determine the properties 
of this critical point, we apply I’Hopital’s Rule a third time to obtain: 

(5.4.10) lim |^(s,x) h^ + l^alci + 120 ^ -F 3 ( 6 ^ - af)^)^^ ^ 

2=1 

Since ^ is uniformly bounded from below and r < 2m -I- 1, from the proof of 


Lemma 4.2.6 it follows that this term can be made negative if |x| is large enough. 
By the point-wise second derivative rule, s = 0 is a local maximum of the function 
(/)(s,x), and therefore, we can find he such that: for any h < he, 4>(h,x) is bounded 
by (j){h, x) < 0, as required by the estimate in (5.4.3). Following the proof 


of Lemma |4.2.6[ we can also prove that the higher order terms in a Maclaurin 
series expansion of 4>(h, x) in powers of h are dominated by a negative term if |x| is 
large enough. Thus, the bound lim/i_j.o <()(h, x) holds as the distance from the origin 
increases. □ 


Remark 5.4.1. We will henceforth assume that r in (5.1.2) satisfies: 

r = m -|- 1 

This hypothesis ensures that the semigroup associated to the approximation is 


Feller (by Lemma 5.2.2), and that the approximation has a stochastic Lyapunov 
function (by Lemma 5.4.1). 
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Next we apply this stochastic Lyapunov function to study a martingale problem 
associated to the Markov process induced by Qq. 


Lemma 5.4.2. Suppose Assumption 4-1 holds on the drift field, let V{x) be the 
Lyapunov function from Lemma \5.4-l\ and let Lp : K" —>■ K. Consider the local 
martingale: 

M^if) = ip{X{t)) - v?(X(0)) - f QqipiX{s))ds, t>0. 

JO 

If ip satisfies the growth condition: 


(5.4.11) 


Qqip^{x) - 2(p{x)Qq(p{x) <V{x) + C 


for some positive constant C and for all x G M", then M‘^{t) is a global martingale. 

The proof of Lemma 4.3.7| can be used to prove Lemma |5.4.2[ It follows from 
Lemma 5.4.2 that Dynkin’s formula (4.1.10) holds for any function satisfying the 
growth condition (5.4.11). Lemma 4.3.8 (with Q replaced by Qq) provides a prac¬ 
tical way to test if a function satisfies this condition. 

5.5. Accuracy 

Here we consider finite-time and long-time accuracy of the approximation. 


Theorem 5.5.1 (Finite-Time Accuracy). Suppose Assumption 4-1 holds. For 
every ip G C^(U"), there exists a positive constant Cq,^t such that 

(5.5.1) ||P,V-^t<7>lloo 

for any t > Q. 


As in the proof of Theorem 4.5.1 the main issue in the subsequent proof is how 
to deal with the unbounded property of the drift field. 

Proof. Defined the global error as: 

(5.5.2) e{x, t) = Ptip{x) — Pl''ip{x) for any cc S M" and t> 0 . 

This global error e(x, t) is a mild solution of the following inhomogeneous linear 
differential equation on the Banach space H;,(M"'): 


(5.5.3) 


e(x, t) — Qqe{x, t) = {L — Qq)Ptip{x) for any x € K" and t > 0 


with initial condition: 

e(x, 0) = 0 for all x e K" . 

Definition |4.5.1| implies this solution can be written in terms of a variation of con¬ 
stants formula: 

(5.5.4) 


z{x,t) = f Pt_,[L - Qq)Ps(p{x)ds . 
Jo 


The desired estimate follows from the proof of Theorem |4.5.1| with Lemma |4.4.1| 
replaced with Lemma 5.3. 1| □ 


Lemma 5.5.2. Suppose Assumption \4-l\ holds on the drift field. Then the 
Markov process generated by Qq has at least one invariant measure. 

Proof. This result follows by the classical Krylov-Bogoliubov argument [nu, 
since the semigroup associated to Qq is Feller and the process generated by Qq has 
a stochastic Lyapunov function. □ 
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Theorem 5.5.3 (Stationary Distribution Accuracy). Suppose Assumption 4-1 
holds. For every ip G C^(R"), we have 

(5.5.5) \U'^{p)-Ilip)\< h^ 

where C^p is defined as: 


c^p = sup n" 


0<h<hc 


Cp pds 


and Cp^p{x) is given in (5.3.11 with f = PsP- 

Proof. In order to obtain the desired estimate, we average the global error 
with respect to the invariant probability measure of the approximation 11^ (which 
exists by Lemma 5.5.2): 

n^(e)= ^ u\z)e{t,z) 


z^Sh 


y (L - Q)Pspds 


= n'’ 


Passing to the limit t —>■ oo, and using ergodicity of the SDE solution, gives the 
desired estimate. 

□ 






CHAPTER 6 


Tridiagonal Case 


For scalar SDEs, the operators given in ^2.3 have an infinite tridiagonal matrix 
representation. We use this representation in this chapter to show that the invariant 
density, mean first passage time, and exit probability of the approximation satisfy 
a tridiagonal system of equations that can be recursively solved. Explicit formulas 
are given for these solutions, which were used to benchmark the simulations of the 


scalar SDEs appearing in ^3.2 i3.3 and ^3.6 Most of these results are classical, 
but their application to numerical methods for SDEs seems somewhat new. 

6.1. Invariant Density 

Let S = {xi}i^z be a collection of grid points in M ordered by index so that: 

. . . < Xi-i <Xi < Xi+i < ... 


The one-dimensional generators in (2.3.6), (2.3.9), and (2.3.15) can all be written 
as a tridiagonal operator Q whose action on a grid function / : S' —>■ K is given by: 

( 6 . 1 . 1 ) {Qf)i = Qi,i+lfi+l ~ 

where we used the shorthand notation: fi = f{xi) and Qi^i+i = Q{xi, Xi+i) for 
grid points Xi,Xi+i G S. Let Q* be the transpose of Q, so that {Q*)i,j = Qj,i for 
all i,j € Z. Recall, that an invariant (not necessarily integrable) density of Q is a 
grid function : S —)■ K that solves: 

(6.1.2) (t^ Z^cz)i — t^i-t-1,2 (^£z)2-t-l (Qi,i-t-l T — (^d)i T ^2 — l,i(^d)i—l — 6 


for every i G Z. A generator of the form (6.1.1) is always V^i-symmetric’ as the 
next Proposition states. 


Proposition 6.1.1. Suppose Vd satisfies (6.1.2) with seed values (Pd)o = 1 cind 


(Pd)i = Qo,i/Qi,o- Then Q is Vd-symmetric, in the sense that 
(6.1.3) foT all i GZ . 

The converse is also true. 


Proof. The z^^-symmetry condition (6.1.3) immediately implies the invariant 


density condition (6.1.2). For the converse, use induction: the hypothesis on the 


seed values implies that (6.1.3) holds for i = 0. Assume (6.1.3) holds for f = j — 1 > 
0. Since Vd satisfies ([6 .1.2[), it follows that 


The inductive hypothesis implies that the right-hand-side of this equation is zero. 


and therefore, (6.1.3) holds for all z > 0. A similar proof works for z < 0 and is 
therefore omitted. □ 
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6. TRIDIAGONAL CASE 


In addition to ^'d-symmetry, a closed form solution to the recurrence relation 


in (6.1.2) is available. 


Proposition 6.1.2 (Formula for Invariant Density). A formula for the solution 

’id (; 

Qk^k+l 


to (6.1.2) with seed values (Pd)o = 1 o-nd (Pd)i = Qo,i/Qi,o is given by 

rid-i 

ifi>0 


(6.1.4) 


(i^d)i = < 


n 

|z|-l 

n 

k=0 


Qk-\-l^k 

Q — k, — k—l 
Q — k—l, — k 


if i < 0 


ifi = 0 


This expression for in terms of the entries of Q can be easily implemented 
in e.g. MATLAB. 


Proof. The proof shows that Vd in ( |6.1.4 ) satisfies p^-symmetry, which by 
Prop. 6.1.1 implies that i/d solves (6.1.2). Evaluating (6.1.2) at i = 0 implies that 


the formula (6.1.4) holds for i = ±1, since (Pd)±i = Qo,±i/Q±i.o- For i > I, 


i’^d)i+i 




Qk, 


k-\-l 


\k^0 


Qk-\-l,k 


\Q 


2 + 1,2 — 


^ 2-1 

n 

Vfc =0 


Qk, 


k+1 


Qk+l^k 


I Qi,i+1 — {k'd')iQi^i^l 


(vd)i 


which implies the desired result for i > 0. A similar result holds for i < —1. □ 

A Markov process on a countable state space is ergodic (or positive recurrent) 
if it is: 

• irreducible; 

• recurrent (or persistent); and, 

• has an invariant probability distribution. 

Thus, we have the following result HO]. 

Theorem 6.1.3 (Sufficient Conditions for Ergodicity). Let X{t) he the Markov 
process with tridiagonal generator Q given by (6.1.1), Markov semigroup {Pij{t)}t>Q 
defined as: 

P,j{t) =¥{X{t) = xj I A(0) = x,) . 


and invariant density Vd whose formula is given in (6.1.4). This process is ergodic 
if it satisfies the following conditions: 

(Cl): For every integer i 

Qi,i-\-i ^ 0 and Qi^i—i ^ 0 
(C2): The following series converges 

Z = ^(Pd)j < oo 

iGZ 
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Let TTd = va,IZ. Then we have that: 

\ Pij{t) — {T:'d)j \ “^0 as t ^ oo for any i G 'Z . 


Here is a simple sufficient condition for (C2) to hold. 


Lemma 6.1.4. Assume that Vd satisfies (6.1.2) with seed values {vd)o = 1 cind 


(ud)i = Qo,i/Qi,o- If the following inequalities are satisfied: 

Qk,k — 1 


lim < 1 

fe-s-oo Qk+l,k 

then Vd is normalizable. 


and lim 


fc->.-oo Qk-l,k 


< 1 


Proof. This lemma is an application of the ratio test, which gives a sufficient 
condition to establish convergence of a series. From (6.1.4) the quotient this test 
uses can be written as a ratio of jump rates: 


(ud)fe+l _ Qk,k+1 
{Vd)k Qk+l,k 


Hence, the ratio test for convergence of the series 
that: 


1 . Qk.k+l 
lim —- 

fc->oo Qk + l,k 


< 1, 


and similarly for i < 0. 


{vd)i reduces to confirming 


□ 


6.2. Stationary Density Accuracy 


As an application of Formula (6.1.4), we obtain explicit estimates for the error 
in (sup norm) of the numerical stationary density of the generator Qc introduced 
in (pX^. 


Proposition 6.2.1. Consider the SDE (2.3.2) with Q = . 


that: 


• U{x) is three times differentiable for all ai S M, and 

sign{x)U'{x) —>■ oo as |a:| —)■ 00 . 

• M(x) is differentiable and M{x) > 0, for all a; € K. 


Let v{x) = exp(—t7(x)). Suppose that the spatial step size is eonstant and that 
Xi = ih for all i G Z. Consider the second-order scheme (2.3.9) with off-diagonal 
entries given by: 


{Qc)i 


,i+l 


1 + Mi+i f jjih\ 

U—2- 




1 Mj + Mi_i 

P 2 


exp 



Let Vc be an invariant density of Qc normalized so that V(.{xif) = v{xf). Then, 
(6.2.1) sup \v{xi) — iZcixi)] < Chf^ 

XiGS 


c = 


:^supu(a;)|a:| sup |[/"'(0|exp 

X 4<|a;| 



|x| sup \u'''{f)\ 

{<|x| 


where 
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6. TRIDIAGONAL CASE 


Proof. First, we apply the test given in Lemma [6.1.4| to check that an invari¬ 
ant density of Qc is normalizable. The hypothesis on U'{x) implies that: 

c)k^k+l 


lim 


— = lim exp = 0 

Ik fc->oo \ ^ / 


k^oo {Qc)k+l,. 

and similarly for fc < 0. Thus, Qc has a stationary density. 

Next, we use formula |6.1.4| to quantify the point-wise error in the approximation 
to the density v{x). Let Vc(x) denote the stationary density of Qc- (This density 
will be rescaled in order to obtain Vc{x).) Substituting the rates appearing in ( |2.3.9 ) 
into formula |6.1.4| and simplifying yields: 

14-1 


Vc{Xi) = < 


exp 


exp 


KI-1 

'^Y.^uu + uu_,)\ if *<0 


fc =0 


1 


if j = 0 


We consider two cases depending on the value of Xi- Set Vc = exp([/o)) so that 

VciXo) = v[xq). 

• When Xi > 0, 

Vc{xi) = exp J U'{s)ds -Uo + = exp {-Ui + e^) , 

where we have introduced 


4 = 


U'{s)ds - - 




I _i / pxk+i 7 \ 4_1 

= E / - 2 

k=0Q'^^>‘ Q k=0 


Use integration by parts twice to obtain the following expression for 6t: 


4 = ^ 


^ p^k+l 


(s-Xfc+l)^- 


U"'{s)ds = -^U'"(4+) 


U, 

12 


( 6 . 2 . 2 ) 


where Xj^^i = {xk + Xk+i)/^ and Xk < 4 — Xk+i- Next we apply the 
elementary inequality: 

|e“ - 11 < el“l - 1 < |a|el“l VaSK 


to obtain 

WiXt) - l^ciXi) \ = iy{Xi) 




“P(-I2EC'"'K 


- 1 


fc =0 


<^v{xi)\xi\ sup |U"'(C)|exp f ^Ixil sup \U'''{i)\ 


12 




12 ' 
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• When Xi < 0, we obtain: 

rO 


V, 


{xi) = exp (^J U\s)ds - t/o + C/t ^ = exp {-Ui + e^) , 


and as in the previous case, 




k=0 


where € [xk-i^Xk]- Applying (6.2.2) we obtain, 


W(x^) - Vcix^) \ = u{Xi) 






fe=0 


<^v{xi)\xi\ sup |C/"'(0|exp ( ^|xi| sup \U'''{^)\ 


12 


^<k^| 


12 ' 


^<k^l 


Taking the sup over Xi in the last inequalities in each case gives the desired 
estimate. □ 


6.3. Exit Probability 

Here we demonstrate how to compute the exit probability or committor func¬ 
tion of the approximation induced by Q. As a side remark, the committor function 
can be used to analytically describe the statistical properties of trajectories between 
two metastable states of an SDE Given ‘reactant’ and ‘product’ 

states at the points a and b (resp.) with a < b, the committor function is defined 
as the probability that the process initiated at a point x S [o, b] reacts in the sense 
that it first hits b before it hits a. Assume that the reactant and product states 
are the left and right endpoints of the grid S and suppose that the grid points are 
arranged so that: 

xq = a <■■■ < xn = b 

Let qd = {((?d)i}^o denote the committor function associated to the generator Q 
which satisfies: 

iQqd)^ = 0 if 0 < t < A^ , 

1 1 n I 1 1 

(qdjo = 0 , (qd)N = 1 • 

The following proposition gives a formula for qd in terms of the reaction rates and 
invariant density Vd of the approximation. 


Proposition 6.3.1 (Formula for Exit Probability). The committor function 

s: 

(Pd)lQl,0 


which solves (6.3.1) can be written as: 

i-l 


(6.3.2) 


iqd)x = < 






0 


ifO<i<N 


ifi = 0 


where Z is a constant selected so that {qd)N = 1 
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6. TRIDIAGONAL CASE 


Proof. Let qa solve the recurrence relation in ( |6.3.1 ), which implies that 

(6.3.3) (Qfz)^) — Qi,i—l{{qd)i (for t ^ 1) 

with seed values {qd)o = 0 and {qd)i = 1- Then q^ = qdl{qd)N solves the same 
recurrence relation with the desired boundary conditions: ((id)^ = 0 and {qd)N = 1- 
Hence, this proof derives a solution to (6.3.3). Multiply ( |6.3.3 ) by {vd)i and use 
z/d-symmetry to obtain: 

i^^d^i+lQi+l^iiiqd^i+l (9^)2) — (Qcz)^— l) 

Sum these equations over the index i from 1 to any fc > 1 and use summation by 
parts to obtain: 


(^cz)i+ll5i+l,i ( (^tz)i+l {^qd)i) — ^ '^^ {^d)iQi,i—l{{qd^i (Qcz)^— l) 

i=l i=l 

fc+1 k 

i—2 2—1 

{j^d)k+lQk+l,kiiqd)k+l — {qd)k) = {k'd)lQlfi{{qd)l — iQd)o) ■ 


Combine this result with a telescoping sum to obtain: 


{qd)k+l = {qd)k+l — {qd)o + {qd)o 


= - {qd)j) + {qd)o 

j=0 

k 

- {qd)o) + {qd)o for fc > o . 


3=0 


{^d)j + lQj + l,j 


Substituting the seed values {qd)o = 0 and {qd)i = 1 gives the solution to (6.3.3). 
The desired formula is obtained by setting q^ = qd/{qd)N ot Z = l/{qd)N in 

(pXl. □ 


6.4. Mean First Passage Time 

A similar recurrence relation can be derived for the mean first passage time 
(MFPT) of the approximation from any point in {xi,,xpi\ to x^ or xj^. Define r as 

T = min{t > 0 I X(t) G {xl^xb} , -A(0) = Xi} . 

Let Ud = {{ud)i}iLo be the MFPT, which satishes Ud{xi) = E2,;(r) and, 

{Qud)i = -1 if 0 < i < Af , 

(6.4.1) ' ' 

Ud[XLj = Ud { XR ) = 0 

The following proposition gives a formula for Ud in terms of the reaction rates and 
invariant density Vd of the approximation. 
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Proposition 6.4.1 (Formula for MFPT). The MFPT which solves (6.4.1) can 
be written as: 

(6.4.2) 




0 


ifO<i<N 
if i = N or i = 0 


where Z is a constant selected so that {ud)N = 0. 


Proof. This proof is similar to the proof of Proposition |6.3.1[ except that 
we do not introduce an auxiliary solution. Multiply (6.4.1) by {i'd)i and use Vd- 
symmetry to obtain: 

(Pfl)l+ll5l+l,l((^tZ^(l)l+l ir^d^i) 1 ((^d)l (lltz)l—l) — {r^d^i 

Sum these equations over the index i from 1 to any fc > 1 and use summation by 
parts to obtain: 


i—1 i—1 i—1 

fc + l k k 


{^d)k-\-lQk-\-l,k{{'^d)k-\-l {'^d)k^ — (^rf) 1 Ql,0((^rf) 1 (^(i)o) ^ ^ i ■ 

Combine this result with a telescoping sum to obtain: 

{Ud)k+1 = {Ud)k+1 - iud)o + {ud)o 


= + 1 + iUd)o 

1=0 


_ / N {t^d)lQl,0 _ 


k j 


(^c/)z 


To enforce the boundary condition at xn = b we solve 

N-l , y ^ N-i j 


(i^d)i 




= 0 


for {ud)i to obtain 


s^N-l sr^j 
2-^j—o 2-^i—i 


{Vd)i 


{Ud)l = 


^N-l {>^d)lQl,0 

2^1=0 


^ (^d)j + lQj + l,j 

which gives the desired formula with Z = {ud)i- 


□ 
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6. TRIDIAGONAL CASE 


By construction, the mean holding time of the generator Qe in (2.3.15) at 


every grid point Xi is equal to the mean first passage time of Y{t) to 
conditional on Y (0) = Xi. Here we show that this property implies that the Markov 
process associated to Qe also preserves the mean first passage time of Y{t) to 
{xk,XjY for any Xk,Xj G S with Xk < Y(0) = Xi < Xj. Abstractly speaking, this 
result manifests the additive property of mean first passage times. 


Lemma 6.4.2. Let Qe denote the generator in (2.3.15) that uses exact jump 


rates and transition probabilities. Let v{x) be the exact solution to: 


(6.4.3) Lv{x) = —1 and v{xk) = v{xj) = 0 

where L is the generator of the SDE and Xk,Xj G S with Xk < xj. Then v{x) 
satisfies: 

QeV{Xi) = -1 , 


for all grid points Xi G {xk,Xj). 


Proof. The unique solution to the boundary value problem (6.4.3) is: 

fx py 


(6.4.4) 


px py 

fix) = - / / 

JXh Jxh 


fir) 


'X, Jx^ M{y)u{y) 
where Ki is the following constant: 


drdy + Ki 


dy 


M{y)v{y) 


v(r) 


K, = 


Jx, Jx, M{y)v{y) 


drdy 


L 


dy 


M(y)v{y) 


Apply Qe on v{x) to obtain: 


{QeX)i 


Zj - Zi_l 1 

- Vi+I - 


1 1 


^i—1 


[{Zi+l 


Vi Zi+1 -z, 1 

-1- Vi-1 

tli _1 

- Zi-i){vi-i - Vi) + {Zi - Zi-i){v,+i - -yi-i)] . 


To simplify this expression, recall that the exact mean first passage time to {xi-i, cci+i)'^ 
is given by: 


(6.4.5) 


u{x) 



M{y)u{y) 


drdy + Ki 



dy 

M{y)v{y) 


where Ki is the following constant: 


py 


fir) 


K, = 


.i M{y)iy{y) 


drdy 


fXi.ll 


dy 


r._i M{y)v{y) 
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To be sure, u[x) equals v[x) when Xk = Xi-i and Xj = Xi+i. From (6.4.4) and 
(6.4.5) it follows that: 




1 


1 


^2+1 ^i—1 


{Zi+i - Zi-i) 


r^i ry 


/(r) 


M{y)v(y) 


drdy 


- {Zi - Zi-i) 


nXi^-i py 


v{r) 


IXi-X JXk 


M{y)v{y) 


drdy 


1 

Ui 


1 

Ui 


r^i ry 


-Ui +Ui + 

- 2i i 


w_ 

M{y)v{y) 


drdy 


Xi-i JXk 

fx^+i ry 


^i+1 Zi—I 


'Xi-i JXk 


M{y)v{y) 


drdy 


r^i ry 


-Ui + Ui + 


v{r) 

r._i Jxi_^ M{y)i'{y) 


drdy 


rXi+i py 


/(r) 


._1 Jx,.^ M{y)v{y) 


drdy 


\ 


r^i+i 


dy 


dy 


,_i M(y)v{y) 


V M{jj)v(y) 


= -1 


as required. Thus, the truncation error of Qe with respect to the exact MFPT v{x) 
is zero. □ 
























CHAPTER 7 


Conclusion 


This paper presented methods to simulate SDEs that permit fixing the spa¬ 
tial step size of the approximation. This feature of the approximation stabilized 
the method for both finite and long-time simulations. To construct solvers with 
this capability, we considered spatial discretizations of the infinitesimal generator 
of the SDE, which constitutes a shift in perspective from the standard approach 
based on time-discretization of realizations of the SDE. In order to obtain scalable 
approximations, we required that these spatial discretizations satisfy a realizability 
condition. Intuitively speaking, this condition is a non-negativity requirement on 
the weights used in the finite difference approximation of partial derivatives, and im¬ 
plied that the spatial discretizations generate a Markov jump process. Realizations 
of this process can be produced without time discretization error using the SSA 
algorithm. Roughly speaking, the SSA produces a simulation of a continuous-time 
random walk with a user-defined jump size, jump directions given by the columns 
of the (local) noise coefficient matrix, and a mean holding time that is determined 
by the generator. Since we used a Monte-Carlo method to simulate the approxi¬ 
mation, our approach applies to problems that practically speaking are beyond the 
reach of standard numerical PDE methods. Moreover, our approach did not require 
discretizing any part of the domain of the SDE problem. Thus, we were able to 
construct realizable discretizations for SDE problems with general domains using 
simple upwinded and central difference methods that are not necessarily tailored to 
a grid, though other types of discretizations (like finite volume methods) are also 
possible. These realizable discretizations were tested on a variety of SDE problems 
with additive and multiplicative noise. 

In the numerical tests, we found that the mean holding time of the approxi¬ 
mation in any state adjusts according to the stiffness of the SDE coefficients. In 
particular, the mean holding time was inversely proportional to the magnitude of 
drift. We also found that once a bound on the spatial step size is set, the mean 
holding time at any state is determined by the discretization of the infinitesimal 
generator. As we illustrated in the Cox-Ingersoll-Ross and Lotka-Volterra SDE 
problems, this spatial step size can be chosen adaptively according to some mea¬ 
sure of the local Lipschitz constant of the drift field of the SDE. These examples also 
showed that the method does not produce moves outside the domain of definition 
of the SDE. This property is a consequence of the spatial updates being defined to 
lie in the domain of the SDE. Numerical tests on non-self-adjoint, planar diffusions 
illustrated that the integrator is stable and accurate for SDEs with locally Lipschitz 
continuous drift fields, unknown stationary distributions, and irregular coefficients 
(internal discontinuities or singularities). We also considered a Brownian dynam¬ 
ics simulation of a cluster of 13 particles immersed in an unbounded solvent. The 
solvent induced interactions between particles, which were long-range and modeled 
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7. CONCLUSION 


using multiplicative noise. This example demonstrated that the spatial step size 
can be set according to a characteristic length scale in the SDE problem - namely, 
the Lennard-Jones radius - in order to avoid exploding trajectories. 

The first part of the theory analyzed the properties of a second-order accurate, 
realizable discretization on a gridded state space with fixed spatial step size. The 
main issues in this analysis were that the state space may not be finite-dimensional, 
and the linear operator associated to the approximation may not be bounded with 
respect to bounded functions. To deal with the unbounded nature of the state space 
and the generator, we used probabilistic tools to show that the approximation is 
stable. In particular, we showed that under a weak dissipativity condition on the 
drift field, a sharp stochastic Lyapunov function was inherited by the approxima¬ 
tion. A key result in this analysis was that the infinitesimal drift condition we 
derived is uniform in the spatial step size of the approximation. We then used Har¬ 
ris theorem to prove that the approximation is geometrically ergodic with respect 
to a unique invariant probability measure. The stochastic Lyapunov function of 
the approximation was also used to solve a martingale problem associated to the 
approximation. Using a global semimartingale representation of realizations, we 
analyzed some strong properties of the approximation, including the complexity of 
the approximation as a function of the spatial step size. To estimate the global 
error of the approximation, we used a continuous-time analog of the Talay-Tubaro 
expansion. This analog was derived using tools from the theory of semigroups pro¬ 
duced by linear operators on Banach spaces. Specifically, we derived an expression 
for the global error using a variation of constants formula. An analysis of this for¬ 
mula revealed that the accuracy of the approximation at finite and infinite-time is 
determined by the accuracy of the generator of the approximation. 

Variable step sizes were considered in the second part of the theory. We ana¬ 
lyzed a second-order accurate, realizable discretization that uses randomly rescaled 
spatial steps. The main issue in this analysis turned out to be lack of irreducibility 
in the standard sense, and also lack of a strong Feller property. These issues were 
not surprising since the approximation is a pure jump process with a countable 
state space, even though it may not be confined to a grid. Therefore the transition 
probabilities of the approximation were not absolutely continuous with respect to 
Lebesgue measure, and the noise driving the approximation did not have a regu¬ 
larizing effect. That said, we were still able to derive a sharp stochastic Lyapunov 
function, prove finite-time weak accuracy, prove existence of an invariant proba¬ 
bility measure, and prove accuracy with respect to equilibrium expectations. The 
latter results required that the semigroup of the approximation be Feller. To obtain 
this property, we mollified the generator of the approximation in such a way that 
its jump rates were bounded and the approximation still had a sharp stochastic 
Lyapunov function. 

Finally, we analyzed the one-dimensional case where the space-discrete gener¬ 
ators admit an infinite, tridiagonal matrix representation. By solving tridiagonal 
systems of equations, we verified the n-dimensional theory on gridded state spaces, 
and obtained explicit expressions for the stationary density, mean first passage 
time, and exit probability of the approximation. These formulas were used in the 
numerical tests to validate the approximations. 
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